An Overview of Multiscale Simulations of Materials 
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Multiscale modeling of material properties has emerged as one of the grand challenges in material 
science and engineering. We provide a comprehensive, though not exhaustive, overview of the 
current status of multiscale simulations of materials. We categorize the different approaches in the 
spatial regime into sequential and concurrent, and we discuss in some detail representative methods 
in each category. We classify the multiscale modeling approaches that deal with the temporal scale 
into three different categories, and we discuss representative methods pertaining to the each of these 
categories. Finally, we offer some views on the strength and weakness of the multiscale approaches 
that are discussed, and touch upon some of the challenging multiscale modeling problems that need 
to be addressed in the future. 
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I. INTRODUCTION 



Some of the most fascinating problems in all fields of 
science involve multiple spatial and/or temporal scales: 
processes that occur at a certain scale govern the behav- 
ior of the system across several (usually larger) scales. 
The notion and practice of multiscale modeling can be 
traced back to the beginning of modern science (see, for 
example, the discussion iri). In many problems of ma- 
terials science this notion arises quite naturally: the ul- 
timate microscopic constituents of materials are atoms, 
and the interactions among them at the microscopic level 
(of order nanometers and femtoseconds) determine the 
behavior of the material at the macroscopic scale (of or- 
der centimeters and milliseconds and beyond), the latter 
being the scale of interest for technological applications. 
The idea of performing simulations of materials across 
several characteristic length and time scales has therefore 
obvious appeal as a tool of potentially great impact on 
technological innovation?^'^''*. The advent of ever more 
powerful computers which can handle such simulations 
provides further argument that such an approach can ad- 
dress realistic situations and can be a worthy partner to 
the traditional approaches of theory and experiment. 

In the context of materials simulations, one can distin- 
guish four characteristic length levels: 

(1) The atomic scale 10~^m or a few nanometers) in 
which the electrons are the players, and their quantum- 
mechanical state dictates the interactions among the 
atoms. 

(2) The microscopic scale (~ 10~^m or a few microme- 
ters) where atoms are the players and their interactions 
can be described by classical interatomic potentials (CIP) 
which encapsulate the effects of bonding between them, 
which is mediated by electrons. 

(3) The mesoscopic scale (~ 10~^m or hundreds of mi- 
crometers) where lattice defects such as dislocations, 
grain boundaries, and other microstructural elements are 
the players. Their interactions are usually derived from 
phenomenological theories which encompass the effects 
of interactions between the atoms. 

(4) The macroscopic scale {'^ 10~^m or centimeters and 



beyond) where a constitutive law governs the behavior 
of the physical system, which is viewed as a continuous 
medium. In the macroscale, continuum fields such as 
density, velocity, temperature, displacement and stress 
fields, etc are the players. The constitutive laws are usu- 
ally formulated so that they can capture the effects on 
materials properties from lattice defects and microstruc- 
tural elements. 

Phenomena at each length scale typically have a corre- 
sponding time scale which, in correspondence to the four 
length-scales mentioned above, ranges roughly from fem- 
toseconds to picoseconds, to nanoseconds to milliseconds 
and beyond. 

At each length and time-scale, well-established and ef- 
ficient computational approaches have been developed 
over the years to handle the relevant phenomena. To 
treat electrons explicitly and accurately at the atomic 
scale, methods known as Quantum Monte Carlo (QMC)^ 
and Quantum Chemistry (QC)^ can be employed, which 
are computationally too demanding to handle more than 
a few tens of electrons. Methods based on density func- 
tional theory (DFT) and local density approximation 
(LDA)*"-^ in its various implementations, while less accu- 
rate than QMC or QC methods, can be readily applied to 
systems containing several hundred atoms for static prop- 
erties. Dynamical simulations with DFT methods are 
usually limited to time-scales of a few picoseconds. For 
materials properties that can be modeled reasonably well 
with a small number of atoms (such as bulk crystal prop- 
erties or point defects), the DFT approach can provide 
sufficiently accurate results. Recent progress in linear 
scaling electronic structure methods^ has enabled DFT- 
based calculations to deal with a few thousands atoms 
(corresponding to sizes of a couple of nanometers on a 
side) with adequate accuracy. Finally, the semi-classical 
tight-binding approximation (TBA), although typically 
not as accurate as DFT methods, can extend the reach 
of simulations to a few nanometers in linear size and a 
few nanoseconds in time-scale for the dynamicsifi. 

For material properties at the microscopic scale. 
Molecular Dynamics (MD) and Monte Carlo (MC) sim- 
ulations are usually performed employing CIP which can 
often be derived from DFT calculations^^'^^. Although 
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not as accurate as the DFT and TBA methods, the clas- 
sical simulations are able to provide insight into atomic 
processes involving considerably larger systems, reaching 
up to ~ 10^ atomg^'^. The time-scale that MD simula- 
tions based on CIP can reach is up to a microsecond. 

At the mesoscopic scale, the atomic degrees of freedom 
are not explicitly treated, and only larger scale entities 
are modeled. For example, in what concerns the me- 
chanical behavior of solids, dislocations are the objects 
of interest. In treating dislocations, recent progress has 
been concentrated on the so-called Dislocation Dynamics 
(DD) approachiliS-iSiii which has come to be regarded 
as one of the most important developments in computa- 
tional materials science and engineering in the past two 
decades^-. Such DD models deal with the kinetics of 
dislocations and can study systems with a few tens of 
microns in size and with a maximum strain ^0.5% for a 
strain rate of 10 sec~^ in bcc metalsi^. 

Finally, for the macroscopic scale, finite-element (FE) 
methods^ are routinely used to examine the large- 
scale properties of materials considered as an elastic 
continuum2i. For example, FE methods have been 
brought to bear on problems of strain-gradient plasticity, 
such as geometrically necessary dislocations^. Contin- 
uum Navier-Stoke equations are also often used to study 
fluids. 

The challenge in modern simulations of materials sci- 
ence and engineering is that real materials usually exhibit 
phenomena at one scale that require a very accurate and 
computationally expensive description, and phenomena 
at another scale for which a coarser description is satis- 
factory and in fact necessary to avoid prohibitively large 
computations. Since none of the methods above alone 
would suffice to describe the entire system, the goal be- 
comes to develop models that combine different methods 
specialized at different scales, effectively distributing the 
computational power where it is needed most. It is the 
hope that a multiscale approach is the answer to such a 
quest, and it is by definition an approach that takes ad- 
vantage of the multiple scales present in a material and 
builds a unified description by linking the models at the 
different scales. Fig. 1 illustrates the concept of a unified 
multiscale approach which can reach the length and time 
scale that individual methods, developed to treat a par- 
ticular scale accurately, fail to achieve. At the same time, 
the unified approach can retain the accuracy that the in- 
dividual approaches provide in their respective scales, al- 
lowing, for instance, for very high accuracy in particular 
regions of the systems as required. As effective theo- 
ries, multiscale models are also useful for gaining physi- 
cal insight that might not be apparent from brute force 
computations. Specifically, a multiscale model can be an 
effective way to facilitate the reduction and the analysis 
of data which sometimes can be overwhelming. Overall, 
the goal of multiscale approaches is to predict the per- 
formance and behavior of materials across all relevant 
length and time scales, striving to achieve a balance be- 
tween accuracy, efficiency and realistic description. 



Conceptually, two categories of multiscale simulations 
can be envisioned, sequential and concurrent. The se- 
quential methodology attempts to piece together a hier- 
archy of computational approaches in which large-scale 
models use the coarse-grained representations with infor- 
mation obtained from more detailed, smaller-scale mod- 
els. This sequential modeling approach has proven ef- 
fective in systems where the different scales are weakly 
coupled. The characteristic of the systems that are suited 
for a sequential approach is that the large-scale variations 
decouple from the small-scale physics, or the large-scale 
variations appear homogeneous and quasi-static from the 
small-scale point of view. Sequential approaches have 
also been referred to as serial, implicit or message-passing 
methods. The vast majority of multiscale simulations 
that are actually in use are sequential. Examples of 
such approaches abound in literature, including almost 
all MD simulations whose underlying potentials are de- 
rived from electronic structure theoryS2i24, usually ah ini- 
tio calculationsiiiiS. One frequently mentioned2SiS& ex- 
ample of sequential multiscale simulations is the work 
of Clementi et alm^ who used QC calculations to evalu- 
ate the interaction of several water molecules; from this 
database, an empirical potential was parameterized for 
use in molecular dynamics simulations; the MD simula- 
tions were then used to evaluate the viscosity of water 
from atomic autocorrelation functions; and finally, the 
computed viscosity was employed in computational fluid 
dynamics calculations to predict the tidal circulation in 
Buzzard's Bay of Massachusetts. 

The second category of multiscale simulations con- 
sists of the so-called concurrent, or parallel, or explicit 
approaches. These approaches attempt to link meth- 
ods appropriate at each scale together in a combined 
model where the different scales of the system are con- 
sidered concurrently and communicate with some type 
of hand-shaking procedure. This approach is necessary 
for systems that are inherently multiscale, that is, sys- 
tems whose behavior at each scale depends strongly on 
what happens at the other scales. In contrast to sequen- 
tial approaches, the concurrent simulations are still rel- 
atively new and only a few models have been developed 
to date. In a concurrent simulation, the system is often 
partitioned into domains characterized by different scales 
and physics. The challenge of the concurrent approach 
lies at the coupling between the different regions treated 
by different methods, and a successful multiscale model 
seeks a smooth coupling between these regions. 

In principle, multiscale simulations could be based on a 
hybrid scheme, using elements from both the sequential 
and the concurrent approaches. We will not examine 
this type of approach in any detail, since it involves no 
new concepts other than the successful combination of 
elements underlying the other two types of approaches. 

There already exist a few review papers and special 
editions of articles on multiscale simulation of materials 
in literaturaSi^i^^iSSi^SiSiiS. A mathematic perspective of 
multiscale modeling and computation is also available22i. 
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The present overview does not aim to provide another 
collection of various multiscale techniques, but rather to 
identify the characteristic features and classify multiscale 
simulation approaches into rational categories in relation 
to the problems where they apply. We select a few illus- 
trative examples for each category and try to establish 
connections between these approaches whenever possible. 
Since almost all interesting material behavior and pro- 
cesses are time dependent, we will address both the issue 
of length-scales and the issue of time-scales integration 
in materials modeling. 

The paper is organized as follows: In Section II we dis- 
cuss in detail representative examples of sequential mul- 
tiscale approaches in the spatial regime. In Section III 
we present examples of concurrent multiscale approaches, 
also in the spatial regime . In Section IV we discuss 
representative approaches that extend time-scales in dy- 
namical simulations. Section V contains our comments 
and conclusions for the applicability of the various ap- 
proaches. The examples presented in this overview to 
some extent reflect our own research interests and they 
are by no means exhaustive. Nevertheless, we hope that 
they give a satisfactory cross-section of the current state 
of the field, and they can serve as inspiration for further 
developments in this exciting endeavor. 



II. SEQUENTIAL MULTISCALE APPROACHES 

Two ingredients are required in order to construct a 
successful sequential multiscale model: 

(i) It is necessary to have a priori and complete knowl- 
edge of the fundamental processes at the lowest scale in- 
volved. This knowledge or information can then be used 
for modeling the system at successively coarser scales. 

(ii) It is necessary to have a reliable strategy for en- 
compassing the lower-scale information into the coarser 
scales. This is often accomplished by phenomenological 
theories, which contain a few key parameters (these can 
be functions) , the value of which is determined from the 
information at the lower scale. 

This message-passing approach can be performed in se- 
quence for multiple length scales, as in the example cited 
in the introductionSl. The key attribute of the sequential 
approach is that the simulation at a higher level critically 
depends on the completeness and the correctness of the 
information gathered at the lower level, as well as the ef- 
ficiency and reliability of the model at the coarser level. 

To illustrate this type of approach, we will present two 
examples of sequential multiscale approaches in some de- 
tail. The first example concerns the modeling of dislo- 
cation properties in the context of the Peierls-Nabarro 
(P-N) phenomenological model, where the lower scale in- 
formation is in the form of the so-called generalized stack- 
ing fault energy surface (also referred to as the 7-surface) , 
and the coarse-grained model is a phenomenological con- 
tinuum description. The second example concerns the 
modeling of coherent phase transformations in the con- 



text of the phase-field approach, where the lower scale 
knowledge is in the form of ab initio free energies, and 
the coarse-grained model is again a continuum model. 



A. Peierls-Nabarro model of dislocations 

Dislocations are central to our understanding of me- 
chanical properties of crystalline solids. In particular, 
the creation and motion of dislocations mediate the plas- 
tic response of a crystal to external stress. While con- 
tinuum elasticity theory describes well the long-range 
elastic strain of a dislocation for length scales beyond 
a few lattice spacings, it breaks down in the immedi- 
ate vicinity of the dislocation core. There has been a 
great deal of interest in describing accurately the dis- 
location core structure on an atomic scale because the 
core structure to a large extent dictates the dislocation 
properties^^iS^. So far, direct atomistic simulation of 
dislocation properties based on GIF has not been sat- 
isfactory because the CIP is not always reliable or may 
even not be available for the material of interest, espe- 
cially when the physical system involves several types of 
atoms. On the other hand, ab initio calculations are still 
computationally expensive for the study of dislocation 
core properties, particularly of dislocation mobility. Re- 
cently, a promising approach based on the framework of 
the Peierls-Nabarro (P-N) model has attracted consider- 
able interest for the study of dislocation core structure 

and mobility36,37,38,39,40,41,42,43,44^45,46,47^ ^his approach 

when combined with ab initio calculations for the ener- 
getics, represents a plausible alternative to the direct ab 
initio simulations of dislocation properties. 

The P-N model is an inherently multiscale framework, 
first proposed by Peierls^* and Nabarrc^^ to incorporate 
the details of a discrete dislocation core into an essen- 
tially continuum framework. Consider a solid with an 
edge dislocation in the middle (Fig. I^J: the solid con- 
taining this dislocation is represented by two elastic half- 
spaces joined by atomic level forces across their common 
interface, known as the glide plane (dashed line). The 
goal of the P-N model is to determine the slip distribu- 
tion on the glide plane, which minimizes the total energy. 
The dislocation is characterized by the slip (relative dis- 
placement) distribution 

f(a;) = u(a;,0+) - u(a;,0-) (1) 

which is a measure of the misfit across the glide plane; 
u(a;,0+) and u(a;,0^) are the displacement of the half- 
spaces at position x immediately above and below the 
glide plane. The total energy of the dislocated solid 
includes two contributions: (1) the nonlinear potential 
energy due to the atomistic interaction across the glide 
plane, and (2) the elastic energy stored in the two half- 
spaces associated with the presence of the dislocation. 
Both energies are functionals of the slip distribution f (x). 
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Specifically, the nonlinear misfit energy can be written as 

/C30 
j{{{x))dx, (2) 
-C30 

where 7(f) is the generalized stacking fault energy sur- 
face (the 7-surface) introduced by Vitek^". The non- 
linear interplanar 7-surface can, in general, be deter- 
mined from atomistic calculations. For systems where 
CIP are not available or not reliable (for instance, it 
is exceedingly difficult to derive reliable potentials for 
multi-component alloys), ab initio calculations can be 
performed to obtain the 7-surface. On the other hand, 
the elastic energy of the dislocation can be calculated 
reasonably from elasticity theory: the dislocation may 
be thought of as a continuous distribution of infinitesi- 
mal dislocations whose Burgers vectors integrate to that 
of the original dislocation^. Therefore, the elastic energy 
of the original dislocation is just the sum of the elastic 
energy due to all the infinitesimal dislocations (from the 
superposition principle of linear elasticity theory) , which 
can be written as 



Uelas 



dx / dx' In ■ 



L dfix) d{{x') 



\x — x'\ dx dx' 



(3) 

where /x, and v are the shear modulus and Poisson's ra- 
tio, respectively. L is an inconsequential constant intro- 
duced as a large-distance cutoff for the computation of 
the logarithmic interaction energyS. Note that the Burg- 
ers vector of each infinitesimal dislocation is the local 
gradient in the slip distribution, dt{x)/dx. The gradient 
of {{x) is called dislocation (misfit) density, denoted by 
p{x). Since the P-N model requires that atomistic infor- 
mation (embodied in the 7-surface) be incorporated into 
a coarse-grained continuum framework, it is a sequential 
multiscale strategy. Thus the successful application of 
the method depends on the reliability of both 7-surface 
and the underlying elasticity theory which is the basis for 
the formulation of the phenomenological theory. 

In the current formulation, the total energy is a func- 
tional of misfit distribution f(a;) or equivalently p{x), and 
it is invariant with respect to arbitrary translation of p{x) 
and f{x). In order to regain the lattice discreteness, the 
integration of the 7-energy in Eq. |(2Jl was discretized and 
replaced by a lattice sum in the original P-N formulation, 



U, 



misfit 



00 

E 



7(f.)Ax, 



(4) 



with Xi the reference position and Ax the average spac- 
ing of the atomic rows in the lattice. This procedure, 
however, is inconsistent with evaluation of elastic energy 
[Eq. ©] as a continuous integral. Therefore the total 
energy is not variational. Furthermore in the original P- 
N model, the shape of the solution f(a;) is assumed to be 
invariant during dislocation translation, a problem that 
is also associated with the non-variational formulation of 
the total energy. 



To resolve these problems, a so-called Semidiscrete 
Variational P-N (SVPN) model was recently developed^, 
which allows for the first time the study of narrow dislo- 
cations, a situation that the standard P-N model can not 
handle. Within this approach, the equilibrium structure 
of a dislocation is obtained by minimizing the dislocation 
energy functional 



U disl — U elastic + Umisfit + Ustress + Kb luL, (5) 



where 



Uelastic = 7;Xij[Keip[ 
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misfit 
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(6) 
(7) 



(8) 



.(1) 



with respect to the dislocation misfit density. Here, pi 

(2) (3) 

and are the edge, vertical and screw components 
of the general interplanar misfit density at the i-th nodal 
point, and 73 (fi) is the corresponding three-dimensional 
7-surface. The components of the applied stress interact- 
ing with the p\^\ ■* and pf\ are t^^^ = (T21, t^^^ = 0-22 
and r^^^ = (T23, respectively. K, Kg and Kg are the pre- 
logarithmic elastic energy factors^. The dislocation den- 
sity at the z-th nodal point is pi = {fi — fi-i)/{xi — Xi-i) 
and Xij is the elastic energy kernel^. 

The first term in the energy functional, UeiasUc is now 
discretized in order to be consistent with the discretized 
misfit energy, which makes the total energy functional 
variational. Another modification in this approach is 
that the nonlinear misfit potential in the energy func- 
tional, Umisfit, is a function of all three components of 
the nodal misfit, f(xi). Namely, in addition to the misfit 
along the Burgers vector, lateral and even vertical misfits 
across the glide plane are also included. This allows the 
treatment of straight dislocations of arbitrary orientation 
in arbitrary glide planes. Furthermore, because the mis- 
fit vector i{xi) is allowed to change during the process 
of dislocation translation, the energy barrier (referred to 
as the Peierls barrier) can be significantly lowered com- 
pared to its corresponding value from a rigid transla- 
tion. The response of a dislocation to an applied stress is 
achieved by minimization of the energy functional with 

respect to pi at the given value of the applied stress, tI^\ 
An instability is reached when an optimal solution for pi 
no longer exists, which is manifested numerically by the 
failure of the minimization procedure to converge. The 
Peierls stress is defined as the critical value of the applied 
stress which gives rise to this instability. 

The SVPN model has been applied to study vari- 
ous interesting material problems related to dislocation 
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phenomenaSiaSiai. One study involved the understand- 
ing of hydrogen-enhanced local plasticity (HELP) in Al. 
HELP is regarded as one of three general mechanisms 
responsible for H embrittlement of metalsr-"^. There has 
been overwhelming experimental evidence in support of 
HELP, but a theoretical foundation was lacking. In order 
to gain understanding of the physics behind the HELP 
mechanism, Lu et al. carried out ab initio calculations 
for the 7-surface of Al with H impurities placed at the 
interstitial sites^. The 7-surface for both pure Al and 
the Al-I-H systems is shown in Fig. O Comparing the 
two 7-surfaces, one finds an overall reduction in 7 energy 
in the presence of H, which is attributed to the change of 
atomic bonding across the glide plane, from covalent-like 
to ionic-like^. 

The core properties of four different dislocations, screw 
(0°), 30°, 60° and edge (90°) have been studied using 
the SVPN model combined with the ab initio deter- 
mined 7-surface. It was found that the Peierls stress 
for these dislocations is reduced by more than an order 
of magnitude in the presence of H*^ , which is compatible 
with the experimental findings that support the HELP 
mechanism^. Moreover, in order to address the exper- 
imental observation for H trapping at dislocation cores 
and H-induced slip planarity, the H binding energy to 
the dislocation cores was calculated^. These calcula- 
tions showed that there is strong binding between H and 
dislocation cores, that is, H is attracted (trapped) to dis- 
location cores which lowers the core energies. More im- 
portantly, the binding energy was found to be a function 
of dislocation character, with the edge dislocation having 
the greatest and the screw dislocation having the lowest 
binding energy. For a mixed dislocation, the binding en- 
ergy increases with the amount of edge component of the 
Burgers vector. These results suggest that in the pres- 
ence of H, it costs more energy for an edge dislocation 
to transform into a screw dislocation in order to cross- 
slip, since the edge dislocation has almost twice the bind- 
ing energy of the screw dislocation'*". In the same vein, 
it costs more energy for a mixed dislocation to transfer 
its edge component to a screw component for cross-slip. 
Therefore the cross-slip process is suppressed due to the 
presence of H, and slip is confined to the primary glide 
plane, exhibiting the experimentally observed slip pla- 
narity. 

A similar approach was applied to the study of the 
vacancy lubrication effect on dislocation motion in Al. 
From this analysis, it was shown that the role of vacan- 
cies is crucial in reconciling the results of Peierls stress 
measured from different experimental techniquesiS. Very 
recently, a multiple plane P-N model has been developed 
to study dislocation phenomena involving more than one 
glide planes, such as dislocation constriction and cross- 
slipii. Finally we should point out that the P-N model is 
just one example of more general cohesive surface models 
that are built upon the idea of limiting all constitutive 
nonlinearity to certain privileged interfaces, while the re- 
mainder of the material is treated via more conventional 



continuum theoriesi. The same strategy can also be ap- 
plied to the study of fracture and dislocation nucleation 
from a crack tip^. 

It is interesting to note that the analysis of 7-surface 
can provide a qualitative understanding of even more 
complex mechanical properties of materials. For exam- 
ple. Rice and coworkers^^ formulated powerful criteria 
for the brittle behavior of materials, by extending the 
Peierls analysis to geometries involving cracks. Based on 
this framework, Waghmare et al.^^-^'^ were able to pre- 
dict which alloying elements could improve the ductility 
of MoSi2 by analyzing the ab initio determined 7-surface 
of the alloys, and comparing the changes induced by al- 
loying to key features of the 7-surface versus the changes 
induced to the surface energy 7^. Remarkably, certain 
predictions of this relatively simple theoretical modeling 
were borne out by subsequent experiments^. 

We have devoted some attention to the description of 
the P-N model and its implementation using ab initio 
7-surfaces, because it is an ideal case of a sequential 
multiscale model: it consists of a well motivated phe- 
nomenological framework, within which the set of atom- 
istically derived quantities is well defined and complete 
(in this case the 7-surface). In this sense, it fulfills all 
the requirements for a coherent and complete multiscale 
model. There are no doubt limitations to it, arising from 
the range of validity of the phenomenological theory, but 
within this range there are no other ambiguities in con- 
structing the multiscale model. Perhaps, its successes, 
some of which we presented above, are owed to this com- 
plete character of the model. 



B. Phase-field model of coherent phase 
transformations 

The structure-properties paradigm is one of the prin- 
cipal pillars in materials science. The term "structure" 
here refers to structures at many different scales, includ- 
ing the atomic scale geometry determined by the crys- 
talline arrangement of atoms, the structure of the de- 
fects that exist in a material, and the structure that 
emerges as a result of the organization of these defects 
into what is referred to as microstructure. Among these 
structures, the microstructure on the scale of microme- 
ters is often directly tied to the mechanical properties 
of materials, and has therefore attracted great interest 
both in terms of scientific understanding and practical 
applications^&^Mii^^. 

Recently, a powerful sequential multiscale approach 
has been put forward for modeling the precipitate 
microstructure and its evolution in multicomponent 
alloysS^iS&, materials which appear in many technological 
applications. The approach is based on the continuum 
phase-field model whose driving forces (free energies) are 
obtained from combined ab initio calculations and the 
mixed-space cluster expansion technique. One interest- 
ing application of this approach concerned the study of 
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precipitation of the 6''(Al2Cu) phase in Cu-Al alloys dur- 
ing thermal aging^^. 

In the phase-field multiscale approach, the nature of 
phase transformation as well as the microstructures that 
are produced are described by a set of continuous order- 
parameter fieldsSi^. The temporal microstructure evo- 
lution is obtained from solving field kinetic equations 
that govern the time-dependence of the spatially inho- 
mogeneous order-parameter fields. Within the diffuse- 
interface description, the thermodynamics of a phase 
transformation and the accompanying microstructure 
evolution are modeled by a free energy that is a func- 
tion of the order-parameter field, or phase field. For a 
structural transformation, the total free energy can be 
written as: 



^tot — ^bulk + ^int 



elast 



(9) 



where J-buik is the bulk free energy, Tinter is the interfa- 
cial free energy, and T^iast is the coherency elastic strain 
energy arising from the lattice accommodation along the 
coherent interfaces in a microstructure. For a microstruc- 
ture described by a composition field c and a set of struc- 
tural order-parameters, Tyi, ?7p, the first two terms of 
Eq. are given by 



:Fbuik + J'^nter = j {/ [c(r) , 7?^ (r)] -f 1 1 Vc(r) | ^ (10) 
+ \ /3.j (p)V,r7p(r)V,r,p(r)}dF 



where f{c,rjp) is the local free energy density^^ and a 
and Pij (p) are the gradient energy coefficients which con- 
trol the width of the diffuse interface. The elastic strain 
energy is obtained from elasticity theory using the homo- 
geneous modulus approximation's. With the total free 
energy of an inhomogeneous system written as a func- 
tion of order-parameter fields, the temporal evolution of 
microstructures during a phase transformation can be ob- 
tained by solving the coupled Cahn-Hilliard equation for 
a conserved field c, and the time-dependent Ginzburg- 
Landau equation for a non-conserved field r]jjHi2^: 



drjp 
dt 



aV c 

oc 



(11) 



(12) 



where M is related to atom mobility and Lp is the relax- 
ation constant associated with the order-parameter rjp. 
As the above equations illustrate, the continuum phase- 
field methodology depends on three input energies: (1) 
bulk free energies of solid solution and precipitate phases, 
(2) precipitate-matrix interfacial free energies, and (3) 
precipitate/matrix lattice elastic energies. Experimen- 
tal determination of these quantities can be difficult and 
problematic. Therefore a physically motivated method 



for accurately determining these quantities is of critical 
importance to predict the microstructure evolution of in- 
terest. In particular, if the quantities can be determined 
from ab initio calculations, the goal of an "a6 initid^ 
modeling of alloy microstructure evolution would be, to 
a great extent, achieved2Si2^. 

Since direct ab initio calculations of free energies are 
either impractical or impossible with currently available 
computational power, a useful method has been devel- 
oped to extend the ab initio energetics to thermodynamic 
properties of alloy systems with hundreds of thousands 
of atoms'^, referred to as the mixed-space cluster ex- 
pansion (CE). In this scheme, energetics from ah initio 
calculations for a number of small unit cell (^ 10 atoms) 
structures are mapped onto a generalized Ising-like model 
for a particular lattice type, involving 2-, 3-, and 4- 
body interactions plus coherency strain energies^S. The 
Hamiltonian can be incorporated into mixed-space Monte 
Carlo simulations of ~ 10^ atoms, effectively allow- 
ing one to explore the complexity of 2^ configurational 
space. As demonstrated by Vaithyanathan et ai, the 
bulk free energy can be obtained from Monte Carlo sim- 
ulations coupled with thermodynamic integration tech- 
niques. The precipitate/matrix interfacial free energies 
can be determined from similar Monte Carlo simulations 
or from low temperature expansion techniques. The elas- 
tic strain energies are of precisely the same form as the 
coherency strain energy used to generate the mixed-space 
CE. Hence, from a combination of ab initio calculations, a 
mixed-space CE approach, and Monte Carlo simulations, 
one can obtain all the driving forces needed as input to 
the continuum phase-field model. The incorporation of 
these energetic properties, obtained from atomistics, into 
a continuum microstructural model represents a bridge 
between these two length scales, and opens the path to- 
ward predictive modeling of microstructures and their 
evolution. 

To illustrate the use of the method, we mention briefly 
the work of Vaithyanathan et al. who studied the prob- 
lem of precipitation of the 9' (AI2CU) phase in Cu-Al al- 
loys. The free energy of the 9' phase is obtained from ab 
initio calculations of the energy at T = K, coupled with 
the calculated vibrational entropy of this phase. The bulk 
free energies of matrix and precipitate phases are then fit 
to the local free energy as a function of order-parameter 
fields in the phase-field model. T = K interfacial ener- 
gies are determined from supercell calculations, both for 
the coherent interface and for the incoherent interface. 
The anisotropy of these interfacial energies is large and 
has been incorporated in the phase-field model. Elastic 
energy calculations for the coherent strain of AI/AI2CU 
{9') and the calculated lattice parameters of each phase 
determine the elastic driving force in this system. Hav- 
ing determined all the necessary thermodynamic input, 
Vaithyanathan et al. were able, for the first time, to 
clarify the physical contributions responsible for the ob- 
served morphology of 9' precipitate microstructure. The 
agreement between the calculated and experimentally ob- 
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served microstructure of 0' in the Al-Cu alloys was ex- 
cellent, confirming the validity of the approach. 

Although the phase-field model is able to predict com- 
plex microstructure evolution during phase transforma- 
tions, it requires as input phenomenological thermo- 
dynamic and kinetic parameters. For binary systems, 
ah initio calculations can provide these parameters for 
the phase-field model, but it is unrealistic to assume 
that such calculations can be used to determine all the 
thermodynamic information for systems beyond ternary. 
Therefore semi-empirical methods, such as CALPHAD 
(calculated phase-diagram) will remain a useful tool in 
such an endeavori^LSiiS. 



C. Other sequential approaches 

Kinetic Monte Carlo (KMC) simulations, coupled with 
atomistically determined kinetic energy barriers, repre- 
sent a powerful class of sequential multiscale approaches. 
For example, a large body of research has been car- 
ried out for surface growth phenomena with KMC sim- 
ulations whose kinetic energy barrier parameters for 
relevant elemental processes are supplied by ah initio 
calculations^SiSi. In an altogether different field, Cai 
et al. have used KMC method to study dislocation 
motion in Si based on the well-established double-kink 
mechanismjS,. In their approach, the dislocation is rep- 
resented by a connected set of straight line segments 
which move as the cumulative effect of a large number 
of kink nucleation and migration processes. The rate of 
these processes is calculated from transition state theory 
with the transition energy barrier having contributions 
from both atomistically determined energetics (double- 
kink formation and migration energy) and elastic inter- 
actions with other dislocation segments as well as from 
the externally applied stress. 

An example of a multiscale approach, in which KMC 
is a key component, employs the so-called level-set 
method^^'*'* for the largest (macroscopic) scale. This ap- 
proach is particularly well suited for the study of epitaxial 
growth, a subject of great importance in microelectronics 
and optoelectronics applications. In the level-set method, 
growth is described by creation and subsequent motion of 
island boundaries. The model treats the growing film as 
a continuum in the lateral direction, but retains atomistic 
discreteness in the growth direction. In the lateral direc- 
tion, continuum equations representing the field variables 
can be coupled to growth through island evolution, by 
solving the appropriate boundary-value problem for the 
field and using local values of this field to determine the 
velocity of the island boundaries. The central idea be- 
hind the level-set method^ is that any boundary curve 
r, such as a step or the boundary of an island, can be 
represented as the set of values if — Q (the level-set) of a 
smooth function ip. For a given boundary velocity v, the 



equation for ip is 

^+vV</. = (13) 

Growth is naturally described by the smooth evolution of 
Lp determined by this differential equation. In the case of 
multilayer growth, the boundaries Tk{t) of the islands are 
defined as the set of spatial points x for which iy9(x, t) = k 
for fc = 0, 1, 2, . . .. The evolution of the level-set function 
If can be obtained by numerically solving Eq. 1)1 3|l using 
non-oscillatory methods^. The key parameters enter- 
ing the model are diffusion constants (the terrace and 
island-edge diffusion constants) which can in principle 
be supplied from atomistic calculations, through the fol- 
lowing procedure (see Fig. first, the atomistic pro- 
cesses are identified which are responsible for terrace or 
island-edge diffusion and their energetics are analyzed us- 
ing atomistic (possibly ah initio) calculations; next, the 
energy barriers for the atomistic processes are incorpo- 
rated in a KMC model which provides the means for 
coarse-graining the atomistic degrees of freedom to a few 
mesoscopic degrees of freedom describing the evolution 
of surface features (the island step edges) ; finally, the re- 
sults of the KMC model are coarse-grained to provide 
the input to the level-set equations, that is, they define 
the values of the boundary velocity v which depends on 
the local surface morphology. The coarse-graining be- 
tween scales eliminates degrees of freedom that are not 
essential, making the passage to the next scale feasible. 
For example, in the illustration shown in Fig. ^ the 
smallest step width in the KMC scale corresponds to a 
two-atom wide region at the microscopic scale, a situ- 
ation that is relevant to the Si(lOO) surface and possi- 
bly to other semiconductor surfaces (such as III-V com- 
pound surfaces). In these cases, surface atoms tend to be 
bound to dimer pairs, which is the essential unit that de- 
termines the step structure, even though the underlying 
dynamics may be determined by the motion of individ- 
ual atoms. Thus, the KMC simulation need only take 
into account structures consisting of dimer units, the dy- 
namics of which determine the step-edge motion needed 
for the level-set simulation. The middle terrace in Fig. 
EJb) is shown as a grid of squares, each representing a 
four-atom cluster and being the minimal unit relevant to 
step motion at the KMC scale in this example, assuming 
that only steps of width equal to two atoms in each di- 
rection are stable. The level-set method is a manifestly 
multiscale approach, combining information from three 
different regimes (atomistic, mesoscopic and continuum) 
into a neatly integrated scheme. Recently, the level-set 
method has also been applied to study dislocation dy- 
namics in alloysSi. 

Yet another sequential multiscale approach has been 
successfully applied to the study of crystal plasticity. 
This is the DD method mentioned earlier, incorporating 
dislocation motion at the macroscopic scale, the mech- 
anism ultimately responsible for crystal plasticity. In 
order to predict the mechanical properties of materials 
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using DD simulations, a connection between micro-to- 
meso scales must be established because dislocation in- 
teractions at close range (when the cores intersect, for 
instance), are totally beyond the reach of continuum 
models. Along these lines, Bulatov et al. were able to 
study dislocation reactions and plasticity in fee metalsSS 
that compare well with deformation experiments, by inte- 
grating the local rules derived from atomistic simulation 
of dislocation core interactions into the DD simulations. 
The same idea has been further explored by Rhee et al. 
in a study of the stage I stress-strain behavior of bcc 
single crystalsi^. 



III. CONCURRENT MULTISCALE 
APPROACHES 

Broadly speaking, a concurrent multiscale approach is 
more general in scope than its sequential counterpart be- 
cause the concurrent approach does not rely on any as- 
sumptions (in the form of a particular coarse-graining 
model) pertaining to a particular physical problem. As 
a consequence, a successful concurrent approach can be 
used to study many different problems. For example, dis- 
location core properties, grain boundary structure and 
crack propagation could all be modeled individually or 
collectively by the same concurrent approach, as long 
as it incorporates all the relevant features at each level. 
What is probably most appealing, however, is that a con- 
current approach does not require a priori knowledge of 
the physical quantities or processes of interest. Thus, 
concurrent approaches are particularly useful to explore 
problems about which little is known at the atomistic 
level and its connection to larger scales, and to discover 
new phenomena. We discuss below three instances of 
concurrent approaches in some detail, and mention some 
additional examples more briefly. 



A. Macroscopic Atomistic Ah initio Dynamics 

Fracture dynamics is one of the most challenging prob- 
lems in materials science and solid mechanics. Despite 
nearly a century of study, several important issues remain 
unsolved. In particular, there is little fundamental under- 
standing of the brittle to ductile transition as a function 
of temperature in materials; there is still no definitive ex- 
planation of how fracture stress is transmitted through 
plastic zones at crack tips; and there is no complete un- 
derstanding of the disagreement between theory and ex- 
periment regarding the limiting speed of crack propaga- 
tion. These difficulties stem from the fact that fracture 
phenomena arc governed by processes occurring over a 
wide range of length scales that are all connected, and all 
contribute to the total fracture energy^. In particular, 
the physics on different length scales interacts dynami- 
cally, therefore a sequential coupling scheme would not 
be adequate for the study of fracture dynamics. 



To address these challenges, Abraham, Broughton, 
Bernstein and Kaxiras developed a concurrent multi- 
scale modeling approach that dynamically couples dif- 
ferent length scaleaSi2&. This multiscale methodology 
aims at linking the length scales ranging from the atomic 
scale, treated with a quantum-mechanical tight-binding 
approximation method, through the microscale, treated 
via the classical molecular dynamics method, and fi- 
nally to the mesoscale/macroscale treated via the finite 
element method in the context of continuum elastic- 
ity. These authors applied this unified approach, termed 
macroscopic-atomistic- a6 initio dynamics (MAAD), to 
the study of the dynamical fracture process in Si, a typ- 
ical brittle material. In traditional studies of fracture, 
only the continuum mechanics level (employing, for in- 
stance, the FE method) is usually invoked to account for 
the macroscopic behavior. But since there is no natu- 
ral small-length cutoff present in the continuum mechan- 
ics approach, any important aspect of the atomic-scale 
mechanisms for fracture is completely missed. This can 
be remedied by introducing classical MD to the simula- 
tions. In particular, the MAAD approach employed the 
Stillinger- Weber-- interatomic empirical potential for Si 
to perform MD calculations at the atomistic level, for a 
large region of the material near a crack tip. However, the 
treatment of formation and breaking of covalent bonds 
at the atomic scale is not reliable with any empirical po- 
tential, because bonds between atoms are an essentially 
quantum mechanical phenomenon arising from the shar- 
ing of valence electrons. On the other hand, small devi- 
ations from ideal bonding arrangements can be captured 
accurately by empirical potentials, because they are to 
first approximation harmonic, a feature that is easily in- 
corporated in empirical descriptions of the interaction 
between atoms. Therefore, it was deemed necessary to 
include a quantum mechanical approach into the simula- 
tions for a small region in the immediate neighborhood 
of the crack tip, where bond breaking is prevalent dur- 
ing fracture, while further away from this region the em- 
pirical potential description is adequate. The particular 
methodology chosen to model the immediate neighbor- 
hood of the crack-tip, a semi-empirical nonorthogonal 
tight-binding schemeSi, describes well the bulk, amor- 
phous, and surfaces properties of Si. Fig. |S1 shows the 
spatial decomposition of the computational cell into five 
different dynamic regions of the simulation: the contin- 
uum FE region at the far-field where the atomic dis- 
placements and strain gradients are small; the atomistic 
MD region around the crack with large strain gradients 
but with no bond breaking; the quantum mechanical re- 
gion (labelled TB because of the use of the tight-binding 
method) right at the crack tip where atomic bonds are 
being broken and formed; the FE/MD hand-shaking re- 
gion; and the MD/TB hand-shaking region. The total 
Hamiltonian, Htot for the entire system was written as: 

Htot = HFE{{viM&FE) + 
HMDi{r.r}eMD) + 
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HTB{{r,r} eTB) + 
HFE/MD{{u,u,r,r} e FE/MD) + 
HMD/TB{{r,r}eMD/TB) 

The degrees of freedom of the Hamiltonian are atomic po- 
sitions r and velocities f for the TB and MD regions, and 
displacements u and their time rates of change u for the 
FE regions. Equations of motion for all the relevant vari- 
ables in the system are obtained by taking appropriate 
derivatives of this Hamiltonian. All variables can then 
be updated in lock-step as a function of time using the 
same integrator. Thus the entire time history of the sys- 
tem may be obtained numerically given an appropriate 
set of initial conditions. Following trajectories dictated 
by this Hamiltonian leads to evolution of the system with 
conserved total energy, which ensures numerical stability. 

The individual approaches at each level (FE, MD and 
TB) are well established and tested methods. What was 
much more important in this study was the seamless 
hand-shaking of the different methods at the interface 
of the respective domains, namely the hand-shaking al- 
gorithms between FE and MD regions and between the 
MD and TB regions. We present here the main ideas 
behind the coupling of the different regions. 
FE/MD coupling: To achieve the FE/MD hand-shaking, 
the FE mesh spacing is scaled down to atomic dimen- 
sions at the interface of the two regions. In Fig. |S1 
the FE nodes are indicated as small open circles con- 
nected by thin lines. Moving away from the FE/MD 
region and deep into the continuum, one can expand the 
mesh size. In this way, the atomistic simulation is em- 
bedded in a large continuum solid, indicated by a green- 
colored region in Fig. ISfa). FE cells contributing fully 
to the overall Hamiltonian (unit weight) are marked with 
thin solid lines, while cells contributing to the hand-shake 
Hamiltonian (half weight) are represented by thin dashed 
lines. Interactions between the atoms on the MD side, 
which are represented by an interatomic potential, carry 
full weight when fully inside the MD region (thick solid 
lines joining neighboring atoms) and half weight (thick 
dashed lines) when they cross the boundary, with one of 
the neighbors effectively represented by a node in the FE 
region. The FE/MD interface is chosen to be far from 
the fracture region. Hence, the atoms of the MD region 
and the displacements of the FE lattice can be unam- 
biguously mapped onto one another. The assignment of 
weights equal to unity within each region and equal to 
one half at the interface is arbitrary and can be general- 
ized by the introduction of a smooth step function. 
MD/TB coupling: At this interface, the atoms treated 
quantum mechanically are shown in red while those 
treated classically are shown in green. The dangling 
bonds at the edge of the TB region are terminated with 
pseudo-hydrogen atoms. The Hamiltonian matrix ele- 
ments of these pseudo-hydrogen atoms are carefully con- 
structed to tie off a single Si bond and to ensure the 
absence of any charge transfer when that atom is placed 
in a position commensurate with the Si lattice. In other 



words, the TB terminating atoms are fictitious mono- 
valent atoms forming covalent bonds with the strength 
and length of bulk Si bonds. These fictitious atoms were 
called "silogens" : they behave mechanically just like Si, 
but chemically like H. The TB Hamiltonian including 
silicon-silicon and silicon-silogen matrix elements is then 
diagonalized to obtain electronic energies and wavefunc- 
tions, from which the total energy can be computed. 
Thus, at the perimeter of the MD/TB region, there are 
silogens sitting directly on top of the atoms of the MD 
region, which are shown as the smaller red circles on top 
of green circles in Fig. |S1 On one side of the TB/MD 
interface, the bonds to an atom are derived from the TB 
Hamiltonian, and are shown as shaded regions in Fig. |S1 
to indicate the electronic distribution responsible for the 
formation of the covalent bonds. On the other side of 
the interface, the bonds are derived from the interatomic 
potential of the MD simulation. The MD atoms of the 
interface have a full complement of neighbors, including 
neighbors whose positions are determined by the dynam- 
ics of atoms in the TB region; these are shown as small 
green circles on top of the red circles in Fig. O As before, 
the TB code updates atomic positions in lock-step with 
its FE and MD counterparts. 

The MAAD approach was employed to study the brit- 
tle fracture of Si in a geometry containing a small crack 
(notch) within an otherwise perfect solid, with the ex- 
posed notch face in the (100) plane and the notch pointed 
in the (010) direction. The system consisted of 258,048 
mesh points in each FE region, 1,032,192 atoms in the 
MD region, and approximately 280 unique atoms in the 
TB region (for computational reasons, the entire region 
modeled by the TB method was broken into smaller, par- 
tially overlapping regions, each assigned to a different 
processor in a parallel implementation). The lengths of 
the MD region are 10.9 A for the slab thickness along 
the front of the crack, 3649 A in the primary direction of 
propagation, and 521 A in the direction of pull (before 
pulling). Periodic boundary conditions were imposed at 
the slab faces normal to the direction of the crack prop- 
agation (along the front of the crack). The wall-clock 
time for a TB force update was 1.5 s, that for the MD 
update was 1.8 s, and that for the FE update was 0.7 s. 
The TB region was relocated after every 10 time steps 
to ensure that it remains at the very tip of the propagat- 
ing crack. The computational slab was initialized at zero 
temperature, and a constant strain rate was imposed on 
the outermost FE boundaries defining the opposing hor- 
izontal faces of the slab. Furthermore, a linear velocity 
gradient was applied within the slab, which results in an 
increasing internal strain with time. It was observed that 
the Si solid failed in brittle fashion at the notch tip when 
the material is stretched by ~ 1.5%. The limiting speed 
of crack propagation was found to be 85% of the Rayleigh 
speed with the chosen computational setup. In the course 
of the simulation, the straight-ahead brittle cleavage of 
the Si slab left behind a rough surface, with increasing 
roughening as a function of crack distance. Based on 
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these results, the authors suggested that the roughening 
surface is due to the spawning of dislocations with low 
mobility on the time-scale of the crack motion. 

A general problem associated with domain decomposi- 
tion, as in the MAAD simulations, is the spurious reflec- 
tion of elastic waves (phonons) at the domain bound- 
aries due to the changes in system description across 
the boundaries. For example, such effects have been ob- 
served in the atomistic modeling of dislocation motion2^, 
crack propagationp22iS^i2S*2&, and energetic particle-solid 
collisions^^-^^, all of which involved some domain cou- 
pling scheme. Since the MAAD method involves do- 
main decomposition into the TB, MD and FE regions, 
the quality of coupling between different regions needs to 
be examined. In a subsequent paper, the same authors 
reported that there was no visible reflection of phonons 
at the FE/MD interface, and no obvious discontinuities 
at the MD/TB interface—-. Thus, in this scheme the cou- 
pling between the various domains is indeed performed in 
a seamless manner, closely mimicking the actual behav- 
ior of the physical system under investigation. Overall, 
the MAAD approach represents the state of the art of 
current multiscale simulation strategies. It is a finite- 
temperature, dynamic and parallel algorithm which, at 
least as far as general computational aspects are con- 
cerned, is applicable to any type of material. 

Ongoing efforts are exploring the possibility of apply- 
ing the MAAD strategy to study chemical effects on me- 
chanical properties of metallic alloys, such as impurity 
effects on dislocation motion, crack nucleation and prop- 
agation in various metals. There is an important qual- 
itative difference between such systems and the study 
of brittle fracture of Si mentioned above: the nature of 
bonds in metallic systems is very different from the sim- 
ple covalent bonds in Si. This makes necessary the devel- 
opment of a different way of coupling the quantum me- 
chanical to the classical atomistic region, because it is no 
longer feasible to terminate the bonds at the boundary of 
the quantum region by simply saturating them with fic- 
titious atoms like the silogens. In such endeavors, other 
more efficient and versatile quantum mechanical formu- 
lations are desirable. One candidate is the linear scal- 
ing real-space kinetic energy functional method^^. This 
method approximates the non-interacting kinetic energy 
of DFT as a functional of electron density, and elec- 
tronic wave-functions are thus eliminated from calcula- 
tions, and therefore the method is termed as orbital- free 
density functional theory (OFDFT). As a consequence, 
no diagonalization of the electronic Hamiltonian and no 
sampling of reciprocal space are necessary, making the 
method computationally efficienliiffi. In particular, the 
explicit real-space feature of this approach makes it nat- 
urally suitable for domain coupling within the MAAD 
framework. While efforts to construct a fully functioning 
scheme along these lines are continuing, we believe this 
is a promising method with great potential for applica- 
tions in metallic systems, which are difficult to handle 
with other techniques. 



B. Quasicontinuum model 

One observation from many large-scale atomistic sim- 
ulations is that only a small subset of atomic degrees 
of freedom do anything interesting. The great major- 
ity of the atoms behave in a way that could be de- 
scribed by effective continuum models like elasticity the- 
ory. The computation and storage of the uninteresting 
degrees of freedom - necessary for a fully atomistic cal- 
culation - consume a large proportion of computational 
resources. This observation calls for novel multiscale ap- 
proaches which can reduce the number of degrees of free- 
dom in atomic simulations^Siiifl^. One such approach 
proposed by Tadmor, Ortiz and Phillips is particularly 
promising and has yielded considerable success in many 
applicationsiS^. This concurrent multiscale approach is 
called the quasicontinuum method, which seamlessly cou- 
ples the atomistic and continuum realms. The chief ob- 
jective of the approach is to systematically coarsen the 
atomistic description by the judicious introduction of 
kinematic constraints. These kinematic constraints are 
selected and designed so as to preserve full atomistic res- 
olution where required - for example, in the vicinity of 
lattice defects - and to treat collectively large numbers of 
atoms in regions where the deformation field varies slowly 
on the scale of the lattice. Variants of the quasicontin- 
uum model have been developed and applied in different 

„^^^^^^^„^-^„1 03.104.105.106.107.108.109.110.111.112.113.114 rpj^^ 

essential building blocks of the static quasicontinuum 
model are: (1) the constrained minimization of the atom- 
istic energy of the solid; (2) the use of summation rules to 
compute the effective equilibrium equations; and (3) the 
use of adaptation criteria in order to tailor the computa- 
tional mesh to the structure of the deformation field. An 
extension of the method to finite-temperature has also 
been proposed^^^. 

The quasicontinuum modcl^^^ starts from a conven- 
tional atomistic description, which computes the energy 
of the solid as a function of the atomic positions. The 
configuration space of the solid is then reduced to a sub- 
set of representative atoms. The positions of the remain- 
ing atoms are obtained by piecewise linear interpolations 
of the representative atom coordinates, much in the same 
manner as displacement fields are constructed in the FE 
method. The effective equilibrium equations are then 
obtained by minimizing the potential energy of the solid 
over the reduced configuration space. A direct calcula- 
tion of the total energy in principle requires the evalua- 
tion of sums that are extended over the full collection of 
atoms, namely, 

N 

Etot = Y.^^^ (14) 

i=l 

where N is the total number of atoms in the solid. The 
full sums may be avoided by the introduction of approx- 
imate summation rules. For example, the lattice quadra- 
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ture analog of Eq. H14|l can be written as 

i=l 

where rii is the quadrature weight that signifies how many 
atoms a given representative atom stands for in the de- 
scription of the total energy, and Ei is the energy of i-th 
representative atom. Note that in this case the sum is 
over the representative atoms only. In the quasicon- 
tinuum approach, the FE method serves as the numeri- 
cal tool for determining the displacement fields, while an 
atomistic calculation is used to determine the energy of 
a given displacement field. The positions of the coarse- 
grained atoms are needed because the energy of the rep- 
resentative atoms depends on them. This approach is 
in contrast to standard FE schemes, where the constitu- 
tive law is introduced through a phenomenological model. 
The selection of the representative atoms may be based 
on the local variation of the deformation field. For exam- 
ple, near dislocation cores and on planes undergoing slip, 
full atomistic resolution is attained with adapted mesh- 
ing. Far from defects or other highly stressed regions, the 
density of representative atoms rapidly decreases, and the 
collective motion of very large numbers of atoms is dic- 
tated, without appreciable loss of accuracy, by a small 
number of representative atoms. 

The quasicontinuum method has been applied 
to a variety of problems, including dislocation 
structures^S2iiS4, interactions of cracks with grain 
boundariesiS^, nanoindentationsiiSiiiSiiii, dislocation 
junctionsiS^, atomistic scale fracture processiS^, etc. By 
way of example, Shenoy et al. applied the method to 
study the interaction of dislocations with grain bound- 
aries (GB) in AliS^. In particular, they considered a 
reformulation of the quasicontinuum model that allows 
for the treatment of interfaces, and therefore of polycrys- 
talline solids. As the first test of the model, they com- 
puted the GB energy and atomic structure for various 
symmetric tilt GB's in Au, Al, and Cu using both direct 
atomistic calculations and the model calculations. They 
found excellent agreement between the two sets of cal- 
culations, indicating the reliability of the model for their 
purpose. In the study of Al, they used nanoindentation- 
induced dislocations to probe the interaction between dis- 
locations and GB's. Specifically, they considered a block 
oriented so that the (111) planes are positioned to allow 
for the emergence of dislocations which then travel to the 
S 21(241) GB located at ^ 200 A beneath the surface 
[see Fig. El^a)]. First, the energy minimization is per- 
formed to obtain the equilibrium configuration of the GB, 
then a mesh is constructed accordingly as shown in Fig. 
ISJa). The region that is expected to participate in the 
dislocation-GB interaction is meshed with full atomistic 
resolution, while in the far fields the mesh is coarser. The 
displacement boundary conditions at the indentation sur- 
face are then imposed onto this model structure, and af- 
ter the critical displacement level is reached, dislocations 



are nucleated at the surface. With the EAM potential- 
supplying the atomistic energies in the quasicontinuum 
approach, they observed closely spaced (15 A) Shockley 
partials nucleated at the free surface. As seen from Fig. 
|Hfb), the partials are subsequently absorbed at the GB 
with the creation of a step at the GB and no evidence of 
slip transmission into the adjacent grain is observed. The 
resultant structure can be understood based on the con- 
cept of displacement shift complete latticeii^ associated 
with this symmetric tilt GB. As the load is increased, the 
second pair of Shockley partials is nucleated. These par- 
tials are not immediately absorbed into the GB, but in- 
stead form a pile-up [Fig. 13b)]. The dislocations are not 
absorbed until a much higher load level is reached. Even 
after the second set of partial dislocations is absorbed 
at the GB, there is no evidence of slip transmission into 
the adjacent grain, although the GB becomes much less 
ordered. The authors argued that their results give hints 
for the general mechanism governing the absorption and 
transmission of dislocations at GB's. The same work also 
studied the interaction between a brittle crack and a GB 
and observed stress-induced GB motion (the combina- 
tion of GB sliding and migration). In this example, the 
reduction in the computational effort associated with the 
quasicontinuum thinning of degrees of freedom is signif- 
icant. For example, the number of degrees of freedom 
associated with the mesh of Fig. Ela) is about 10"*, three 
orders of magnitude smaller than what would be required 
by a full atomistic simulation (10^ degrees of freedom). 

Recently, the quasicontinuum model has been ex- 
tended to complex Bravais lattices"'^^'' whereby more com- 
plicated materials can be handledii^. But because of the 
expression for the total energy adopted in Eq. H14|) . 115|l . 
the actual atomistic methods that can be implemented in 
the quasicontinuum model are limited to ones that can 
be easily cast in such a form, if one insists on having 
the ability to resolve the FE nodes all the way to the 
atomic scale. This limit is often referred to in the lit- 
erature as the "non-local" regime of the quasicontinuum 
method. In contrast, the "local" limit refers to the case 
where each FE node represents a very large number of 
atomistic degrees of freedom, which is modeled as cor- 
responding to an infinite solid homogeneously deformed 
according to the local strain at the node. In this limit, 
it is advantageous to employ effective Hamiltonians to 
compute energetics for the representative atoms. Such 
Hamiltonians can be constructed from ab initio calcula- 
tions, and may include physics that atomistic simulations 
based on classical interatomic potentials (such as EAM) 
are not able to capture. For example, by constructing 
an effective Hamiltonian parametrized from ab initio cal- 
culations, Tadmor et al. were able to study polarization 
switching in piezoelectric material PbTiOs in the con- 
text of the quasicontinuum model in the local limifeiiS. 
This particular Hamiltonian includes the following terms: 
the elastic energy of the lattice, the coupling between 
strain and atomic displacement, harmonic and anhar- 
monic phonon energy contributions, the interaction of 
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atomic displacement with the applied electric field, and 
the electrostatic energy. With this effective Hamiltonian, 
it was shown that the behavior of a large-strain ferroelec- 
tric actuator under the application of competing external 
stress and electric fields can be modeled successfully, re- 
producing, for example, all the important features of the 
experimental strain vs. electric field curve for the actu- 
ator. The advantage of these simulations is that they 
also provide insight into the microscopic mechanisms re- 
sponsible for the macroscopic behavior, making possible 
the improvement of design of the technologically useful 
materials. 

One pitfall of the quasicontinuum model is the so- 
called "ghost force" at the interface between the local 
region, identified with slow variation of the deforma- 
tion gradient, and the nonlocal region, identified with 
rapid variation of the deformation gradienliiSS. The er- 
ror arises from the discontinuity between the neighbor- 
ing cells where the cell sizes are less than the range of 
the atomistic potential. Care must be taken to cor- 
rect these "ghost forces"ii2^. Finally we should point out 
that the quasicontinuum approach also shares certain fea- 
tures with sequential approaches, namely, the constitu- 
tive equation for the FE nodes is drawn from atomistic 
calculations (akin to message passing in sequential ap- 
proaches). The reason we categorize it as a concurrent 
multiscale approach is that the atomistic and FE calcula- 
tions are performed concurrently rather than in sequence, 
because the range of deformations encountered in various 
parts of the system are not know beforehand. Moreover, 
some sort of domain partitioning (meshing) is involved 
in the quasicontinuum approach. 



C. Coarse-grained molecular dynamics 

Mesoscopic elastic systems, and in particular micro- 
electro-mechanical-systems (MEMS), recently have cap- 
tured a great deal of attention and research interest as 
micro-machines and devices. However, there is serious 
concern regarding their mechanical integrity and stabil- 
ity in applications because these sub-micron devices are 
so minuscule that structural defects and surface effects 
could have large impact on their performance. On the 
other hand, the computational study of the mechani- 
cal properties of the MEMS has turned out to be ex- 
tremely difficult because they are too small in size for 
finite-element simulations (at the limit where continuum 
elasticity theory may be no longer valid), but too large 
for atomistic simulationsiiiii2£. To resolve this prob- 
lem, a concurrent multiscale simulation strategy called 
coarse-grained molecular dynamics (CGMD) has been 
developed by Rudd and Broughtonii2ii22,. This approach 
bears some resemblance with the quasicontinuum model, 
yet there exist important differences between the two. 

The CGMD approach is based on a statistical coarse- 
graining prescription. In particular, the model aims at 
constructing scale-dependent constitutive equations for 



different regions in a material. In general, the material of 
interest can be partitioned into cells, whose size varies so 
that in important regions a mesh node is assigned to each 
equilibrium atomic position, whereas in other regions the 
cells contain many atoms and the nodes need not coincide 
with atomic sites. The CGMD approach produces equa- 
tions of motion for a mean displacement field defined at 
the nodes by defining a conserved energy functional for 
the coarse-grained system as a constrained ensemble av- 
erage of the atomistic energy under fixed thermodynamic 
conditions. The key point of this effective model is that 
the equations of motion for the nodal (mean) fields are 
not derived from the continuum model, but from the un- 
derlying atomistic model. The nodal fields represent the 
average properties of the corresponding atoms, and equa- 
tions of motion (in this particular case Hamilton's equa- 
tions) are constructed to describe the mean behavior of 
underlying atoms that have been integrated out. 

One important underlying principle of CGMD is that 
the classical ensemble must obey the constraint that the 
position and momenta of atoms are consistent with the 
mean displacement and momentum fields. To be specific, 
let the displacement of atom /i be = - x^o where 
x^o is its equilibrium position. The displacement of mesh 
node j is an average of the atomic displacements 

where /j^ is a weighting function, a microscopic analog 
of the FE interpolating functions. Note that Latin in- 
dices, J, k, denote mesh nodes and Greek indices, fi, v, 
denote atoms. A similar relation holds for the momenta 
p^. Since the nodal displacements are fewer or equal 
to the atomic positions in number, fixing the nodal dis- 
placements and momenta does not necessarily determine 
the atomic positions entirely. Therefore some subspace 
of phase space remains not sampled, which corresponds 
to the degrees of freedom that are missing from the mesh 
points. The coarse-grained energy is defined as the aver- 
age energy of the canonical ensemble on this constrained 
phase space: 

E{Uk,Uk) = {HMD)ukMk 

(17) 

where f3 — \/{kBT) is the inverse temperature and Z 
is the partition function. The 3-D delta functions 5{u) 
enforce the mean field constraint [Eq. (|16|l ] . 

When the mesh nodes and the atomic sites are identi- 
cal, the CGMD equations of motion agree with the atom- 
istic equations of motion. As the mesh size increases 
some short-wavelength degrees of freedom are not sup- 
ported by the coarse mesh. But these degrees of freedom 
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are not neglected entirely, because their thermodynamic 
average effect has been retained. This approximation is 
expected to be good if the system is initially in thermal 
equilibrium, and the missing degrees of freedom only pro- 
duce adiabatic changes to the system. The Hamiltonian 
was derived originally for a monoatomic harmonic solid, 
but can be easily generalized to polyatomic solids^^^. Af- 
ter deriving the equations of motion from the assumed 
Hamiltonian for a particular system, one can perform 
the CGMD for the nodal points. 

As a proof of principle, the method was applied to 
one-dimensional chains of atoms with periodic bound- 
ary conditions where it was shown that the CGMD gives 
better results for the phonon spectrum of the model sys- 
tem compared to two different FE methodsiiS. A variety 
of other calculations have also been performed with the 
CGMD to validate its effectivenessi^S-iSi. 

Although the CGMD has proven to be reliable in the 
description of lattice statics and dynamics, the imple- 
mentation of the model varies from system to system. 
This is because different approximations have to be made 
to the Hamiltonian that represents a particular system. 
On the other hand, such approximations can be esti- 
mated and controlled in the CGMD method. This ad- 
vantage makes the CGMD method a good candidate for 
replacing the FE method in the MAAD approach when 
a high quality of FE/MD coupling is required. As we 
alluded earlier, the CGMD approach resembles the qua- 
sicontinuum model in the sense that both approaches 
adopt an effective field model, an idea rooted in the 
renormalization group theory for critical phenomena. In 
both approaches, less important (long wave-length) de- 
grees of freedom are removed while the effective Hamil- 
tonian is derived from the underlying fine-scale (atom- 
istic) model. Although both approaches are developed 
to couple FE and atomistic models, the quasicontinuum 
method is mainly applicable to zero temperature calcula- 
tions while the CGMD is designed for finite temperature 
dynamics. The quasicontinuum model has shown its suc- 
cess in many applications, but the CGMD approach has 
yet to show its wider applicability and versatility. 



D. Other works 

Recently a more general model for the dynamics of 
coarse-grained multiscale systems was proposed by Cur- 
tarolo and CederiiSSi. The model is similar to the 
Migdal-Kadanoff approach in the renormalization group 
theory where the system is coarse-grained through 
a bond moving process. The new potentials are con- 
structed to assure that the partial partition function of 
the system remains unchanged. The information re- 
moved from the coarse-graining process can be quanti- 
fied by the entropy contribution of each step. Although 
the model is shown to produce excellent results for me- 
chanical and thermodynamical properties compared to 
the non-coarse-grained system, so far it is limited to two 



dimensional systems, and its generalization to three di- 
mensions is yet to be achieved and tested. 

Another interesting approach has been developed by 
Shilkrot, Miller and Curtin aimed at linking an atomistic 
region to a "defected" dislocation dynamics region^^"*. In 
this coupled atomistic and discrete dislocation (CADD) 
method, the fully atomistic region is directly coupled to 
a linear clastic continuum region containing dislocations 
which are modeled as continuum elastic line defects. 
The dislocations at the continuum region are treated 
with the standard discrete dislocation methodi^^, and 
the atomistic region may have any kind of atomic scale 
defects. The key aspect of the CADD method is that 
the dislocations can pass between the atomistic and 
continuum regions smoothly. Two developments have 
been made to achieve this goal: 

(1) detection of the dislocation near the 
atomic/continuum interface; and 

(2) a procedure for moving the "right" dislocations 
across the interface. 

So far, this approach has only been implemented in 2D 
systems, but it has been shown to agree quite well with 
the 2D atomistic calculations for Al. 

Some other concurrent approaches are similar to the 
MAAD method but concentrate on linking two differ- 
ent length scales rather than three. For example, Bern- 
stein and Hess'^ have simulated brittle fracture of Si by 
dynamically coupling empirical-potential MD and semi- 
classical TB approaches. In a similar vein, Lidorikis et 
al. have studied stress distribution in Si/Si3N4 using a 
hybrid MD and FE approachiSi. More recently, a first- 
principles Green's function boundary condition method 
has been developed to self-consistently couple the strain 
field produced by a line defect to the long range elastic 
field of the host latticeJ^. 

Concurrent multiscale ideas have also been applied 
to the modeling of biomolecules. In particular, the 
hybrid quantum mechanical and molecular mechanical 
(QM/MM) methods have been gaining ground in the 
study of proteins and enzymes in which the small part 
of a molecule (active site) is modeled by ab initio meth- 
ods while the rest of the molecule can be dealt with by 
a more approximate classical force field theorjJSi. One 
particular implementationi^ of the QM/MM strategy 
is to combine the quantum mechanical self-consistent- 
charge density-functional-based TB method^^^ with the 
CHARMM molecular force fields^^^. This approach has 
been used to study the reactions catalyzed by triosephos- 
phate isomerase and the dynamics of small peptide he- 
lices in wateriiSS. 

Finally we wish to comment that many concurrent 
models such as the ones discussed above are designed for 
covalently-bonded systems. These methods take advan- 
tage of the localized electron bonding across the domain 
interface (between TB/MD and between QM/MM), and 
partition the bonding energy approximately with a cer- 
tain degree of empiricism. But for metallic systems, the 
bonds are not localized or associated with a distinct pair 
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of atoms, therefore the concept of "bonding energy par- 
tition" across the domain interface becomes invaUd, and 
new concepts are needed. Recently several groups have 
exploited the idea of "embedding potential" in simula- 
tions where a region (I) with more accurate description of 
the physics is embedded into another region (II) with less 
accurate description. The influence of region II on region 
I is described by the "embedding potential" which corre- 
sponds to a local one-electron operator in the framework 
of the DFT"3ii34ii35^i36^ Pqj. example, in an effort to 
improve the LDA/DFT description of molecular adsorp- 
tion on surfaces, a coupling method was developed where 
a more accurate (quantum chemical) region (I) was em- 
bedded in a less accurate LDA/DFT region (ji)mJMJM. 
The "embedding potential" is defined as the functional 
derivative of the coupling energy with respect to the elec- 
tron density p/(r) in region I. The total electron density 
Ptot — Pi + Pi I, where pn is the electron density in the 
LDA/DFT region, can be obtained by just LDA/DFT 
calculations for the entire system since the electron den- 
sity is usually well represented by LDA/DFT. ptot is then 
held fixed during the subsequent calculations. By em- 
ploying the OFDFT methodiSS for the coupling energy, 
the "embedding potential" can be explicitly evaluated for 
any given pi{v). The "embedding potential" as an effec- 
tive local one-electron operator can in turn be added to 
the Hamiltonian of region I, and the new electron den- 
sity p/(r) is thus determined. In this way, pi{v) can be 
updated self-consistently for the given ptot- The same 
"embedding potential" idea can be applied to the cou- 
pling between two different DFT regions, or between two 
regions where one is treated with DFT, and the other is 
treated classically, for instance with EAM. 

This last approach, currently under development in 
our group, deserves some elaboration. This approach 
strives to combine quantum mechanics via OFDFT, clas- 
sical mechanics via EAM, and continuum mechanics via 
the quasicontinuum method in a unified description for 
metallic systems. Since the electron density defined in 
the EAM potential along with the EAM nuclei, could 
generate an "embedded potential" that the OFDFT elec- 
trons experience, the coupling energy between OFDFT 
and EAM regions can be explicitly calculated. Further- 
more, the EAM atomistic region can be easily coupled to 
the continuum region based on the nonlocal description 
of the quasicontinuum framework. 



IV. EXTENDING TIME-SCALES 

A. Accelerated dynamics 

As we have seen, MD plays a critical role in modeling 
of materials problems because MD simulations can fol- 
low the actual dynamical evolution of the system with- 
out assuming any mechanism or pathway for the dynam- 
ics, in contrast to, say, MC or molecular statics simula- 
tions. However, MD is typically limited to a time-scale 



of nanoseconds because standard MD simulations follow 
the individual vibrations of all the atoms, whose vibra- 
tion frequencies are of the order of 10^'* s^^. This is 
particularly troublesome for the complex systems whose 
dynamics is characterized by the occurrence of rare but 
important events, such as chemical reactions, diffusion 
processes, and conformational changes. In these systems, 
the existence of energetic barriers much larger than ksT 
that separate the initial from final state, leads to reac- 
tion times far greater than those that can be currently 
accessed computationally. The other reason for extend- 
ing time-scales is that time is a sequential object, and 
the current progress in parallel computing has little im- 
pact on solving the problem. Therefore, algorithms that 
could address the time-scale problem could revolutionize 
the field of computational materials science and engineer- 
ing. 

In the past few years, significant progress has been 
made to accelerate MD simulations. A class of ac- 
celerated dynamics methods, including hyperdynam- 
ics, parallel replica dynamics, and temperature accel- 
erated dynamics, have been developed by Voter and 
coworkersiSLi^Sii^. Although each method accom- 
plishes this acceleration in a different way, transition 
state theory (TST) provides the common theoretical 
foundation. TST is an elegant theory with extensive 
applications in materials sciencei^Sti^iii^Siiiiiiiiii^. In 
TST theory, a state-to-state transition rate constant 
{K'^^'^) is approximated as the fiux through a dividing 
surface in phase space separating the initial and final 
states. A common and useful approximation to TST is 
the harmonic TST in which one assumes that the po- 
tential energy surface (PES) near the minimum can be 
expanded with the harmonic approximation. Thus the 
TST rate constant (the flux through the saddle plane) 
becomes 

K"^^^ = voeicp{-E/kBT), (18) 

where 

-0 - (19) 

n 

i 

Here E is the activation energy (energy difference be- 
tween the minimum and the saddle point), Vi is the 
ith normal mode frequency at the minimum, and v[ is 
the non-imaginary normal model frequency at the sad- 
dle point The analytic integration over the whole 
phase space yields the well-known Arrhenius tempera- 
ture dependence. It is worthwhile to point out that al- 
though the exponent depends only on the barrier height, 
there is no assumption that the trajectory passes exactly 
through the saddle point. For systems where there is 
no re-crossing of the dividing surface and the modes are 
truly harmonic, the rate equation H18|l is exact. The un- 
derlying concept in the accelerated dynamics methods is 
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that the system trajectory is simulated to find an ap- 
propriate pathway for escape from an energy weU by a 
process which takes place much faster in the simulation 
than it would with direct MD. We provide below an ele- 
mentary description of this concept. 

The general formulation of TST rests on two assump- 
tions in order to treat infrequent events: 

(a) it is known in advance what the different equilibrium 
states of the system will be; and 

(b) it is possible to construct a reasonable dividing sur- 
face along the boundaries between initial and final states 
(or equivalently, all the saddle points can be identified). 
Unfortunately the knowledge of states through which 
a system may evolve in most cases (especially in 
complex systems) is incomplete. The hyperdynamics 
method^^^"^^^ is designed to accelerate MD simulations 
without any advance knowledge of either the location of 
the dividing surface or the states through which the sys- 
tem may evolve. Based on TST, Voter has shown that 
it is possible to modify the PES of the system in such a 
way that a simulation on this modified surface exhibits 
the correct relative probabilities of transitions, but with 
enhanced overall transition rates for the system escaping 
from one equilibrium state to the various nearby equilib- 
rium states. The key of this approach is to construct a 
bias potential to raise the energy of the system in regions 
other than at the dividing surfaces. Dynamical evolution 
with the biased potential leads to accelerated transition 
from one equilibrium state to another equilibrium state, 
while the elapsed time is related to statistical properties 
of the system. More precisely, the total time advance for 
a hyperdynamics simulation after n integration steps is 

TL 

where AImd is the time advance for a regular MD tra- 
jectory, AV{r) is the bias potential, and T is the tem- 
perature. The overall computational speed-up is given 
by the average boost factor {thyper /tMo), divided by the 
extra computational cost of calculating the bias potential 
and its derivatives. The evolution of hyperdynamics from 
state to state is correct because the bias potential does 
not change the relative TST rates for different escape 
paths from a given state. The long-time dynamics of the 
simulations are exact to the extent that the dynamical 
corrections to the TST are negligible. Recently, Voter 
has shown that the bias potential and its derivatives can 
be computed in 0(N) fashion without ever constructing 
the HessianiSS. Thus, the implementation of the hyper- 
dynamics method requires only first derivatives of the 
interatomic potential, as for normal MD simulations. 

To demonstrate the effectiveness of the method. Voter 
has studied the diffusion of a Ag adatom on the Ag (100) 
surface at 400 K using an EAM potential for the energet- 
ics. He found that a 3.7 x 10^ steps of hyperdynamics 
run gave an average boost of 1356, for a total time of 
9.89 ± 0.5 /xsec. Each hyperdynamics step required ^ 



30 times the computational time of a direct MD step, 
therefore the net computational boost was a factor of 
45. The rate constants obtained from the calculations 
are in agreement with the full harmonic rate. For a more 
complex system with a 10-atom Ag cluster on the Ag 
(111) surface at 300 K, he achieved an average boost of 
8310 with a hyperdynamics run for 221.2 /xsec. With 
this approach, one should be able to observe novel dif- 
fusion mechanisms that can not be accessed by normal 
MD simulations. 

In order to take advantage of recent advances in par- 
allel computation. Voter proposed the so-called parallel 
replica dynamics methodic to treat infrequent events. 
For a system in which successive transitions are uncorre- 
lated, running a number of independent MD trajectories 
in parallel gives the exact dynamical evolution between 
the states. For a system with correlated crossing events, 
the state-to-state transition sequence is still correct, but 
care must be taken to eliminate or reduce the error as- 
sociated with the simulation time. The parallel replica 
method represents the simplest and most accurate of the 
accelerated dynamics techniques, with the only assump- 
tion being that of infrequent events obeying first-order 
kinetics. To be more specific, the probability distribution 
for the waiting time before the next escape is assumed to 
be 

p{t) = Kel^^{~Kt), (21) 

where K is the rate constant for finding the next escape 
path from the current state. In a system that exhibits no 
correlated crossing events, K is exactly the TST rate con- 
stant (K^^^). In a more general case, in which correlated 
crossings occur, K < K^^^ . For an iV-atom system in a 
particular equilibrium state (potential energy basin), the 
entire system is replicated on each of M available par- 
allel or distributed processors. After a short de-phasing 
stage during which momenta are periodically random- 
ized to eliminate correlations between replicas, each pro- 
cessor carries out an independent constant-temperature 
MD trajectory for the entire A^-atom system, thus ex- 
ploring the phase space within the particular basin M 
times faster than a single trajectory would. Whenever a 
transition is detected on any processor, all processors are 
alerted to stop. The simulation clock is advanced by the 
accumulated trajectory time summed over all replicas, 
i.e., the total time spent exploring phase space within 
the basin before the escape pathway is found. The paral- 
lel replica method also correctly accounts for correlated 
dynamical events where TST is no longer valid. This 
is accomplished by allowing the trajectory that made 
the transition to continue on its processor for a fur- 
ther amount of time Atcorr during which re-crossings 
or follow-on events may occur. The simulation clock is 
then advanced by Atcorr, the final state is replicated on 
all processors, and the whole process is restarted. This 
overall procedure then gives exact state-to-state dynam- 
ical evolution because the escape times obey the correct 
probability distributioniSi. With this approach, signifi- 
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cant extensions of MD time-scales can be achieved. For 
example, in MD simulations of vacancy diffusion on the 
Cu(lOO) surface at 500 K, a 15-processor parallel com- 
puter can give a 14-fold increase in simulation time per 
wall-clock time. Moreover, the parallel replica method 
can be combined with other accelerated dynamics meth- 
ods, such as hyperdynamics to give a multiplicative effect 
in the MD time-scale gaini^. 

In the temperature-accelerated dynamics (TAD) 
method, the transition from state to state is acceler- 
ated by increasing temperaturei^. The transitions that 
should not have occurred at the original temperature are 
then filtered out. The TAD method is more approximate 
than the previous two methods due to the fact that it 
relies on the harmonic TST approximation, but it often 
gives substantially bigger boost than the hyperdynamics 
or the parallel replica dynamics in systems where the ap- 
proximation is justified. Consistent with the accelerated 
dynamics concept, the trajectory in TAD is allowed to 
wander on its own to find each escape path, so that no 
prior information is required about the nature of the reac- 
tion mechanismsi^. Like hyperdynamics, TAD can also 
be combined with the parallel replica method to achieve 
an even higher acceleration on parallel computers. 

B. Finding transition pathways 

As stated earlier, the problem of finding transition 
pathways for infrequent events between two known equi- 
librium (stable or metastable) states is of considerable 
interest. The accelerated dynamics methods are de- 
signed to find the real dynamic pathways via effective MD 
simulations. Other methods requiring no preconceived 
mechanism or transition state have also been developed 
to locate transition pathways. For example, Elber and 
Karpluai^ developed a method to find the transition 
pathway by minimizing the average value of the poten- 
tial energy along the path rather than trying to find the 
path with the lowest barrier. A more popular approach, 
similar in spirit to the Elber-Karplus method, is the so- 
called nudged elastic band (NEB) approacbi^SiiS, which 
focuses on the global character of the path rather than 
on local properties of the PES. 

The NEB method is based on the "chain-of-states" idea 
where a number of images (or replicas or "states") of 
the system are connected together between the end-point 
configurations to trace out a transition pathway If the 
images are connected with springs of zero natural length, 
one can define the object function for this so-called plain 
elastic band (PEB) method in the following: 

5^^^(i?i, i?p_i) = ^2 ^(^^) + E ^(^» - ^-i)'' 

(22) 

where vector R represents the coordinate of the system, 
V is the potential energy of the system, and k is the 
spring constant. The spring is introduced to ensure that 



the images are evenly spaced along the path. One would 
envision to find the transition pathway by minimizing the 
object function in Eqn. (|22(l with respect to the interme- 
diate images while keeping the end-point images, Rq and 
Rp fixed. The force acting on the image i is 

= ~VV{R^) + f: (23) 

where 

F,' = fc,+i(i?,+i - i?.) - h{R, - 7?,_i). (24) 

However, as demonstrated by Jdnsson et al., the PEB 
method fails to provide the transition pathway in most 
situationsi^^ . For example, if ki is too large, the elastic 
band becomes too stiff, the transition path would then 
"cut the corner" and thus miss the saddle point region. 
On the other hand, if fcj is small, the elastic band comes 
closer to the saddle point, but the images slide down from 
the energy barrier (avoid the saddle point) , therefore re- 
ducing the resolution of the path in the most critical 
region. Furthermore, by noticing the analogy between 
the object function in the continuum limit and the ac- 
tion of a classical particle of unit mass moving on the 
inverted PES, Jonsson et al. argued that the particle 
would move through the saddle point region with a finite 
velocity affected by the force component perpendicular to 
the curved path. In other words, the images would de- 
viate from the minimum energy path (MEP). The prob- 
lem with "corner cutting" is due to the component of 
the spring force that is perpendicular to the path, which 
tends to pull images off the MEP. The problem with slid- 
ing down results from the component of the the potential 
force or of the true force, VV(i?i) that is parallel to the 
path. The distance between images becomes uneven so 
that the net spring force can balance the parallel compo- 
nent of the true potential force. To cure these problems, 
the NEB method projects out the perpendicular compo- 
nent of the spring force and the parallel component of the 
potential force relative to the path. The force on image 
i becomes 

i^°--VV(i?.)U + i^'-r||f||, (25) 

where f|| is the unit tangent to the transition path and 

VV(^,)U = VV(^0 - VV{Ri) ■ -ri|f||. These force pro- 
jections ("nudging") decouple the dynamics of the path 
itself from the discrete representation of the path. The 
spring force thus does not interfere with the relaxation 
of the images perpendicular to the path, and the relaxed 
configuration of the images satisfies WV{Ri)\± = 0, i.e., 
they he on the MEP. 

The implementation of the NEB in a MD program 
is quite simple. First, the energy and its gradient 
are evaluated for each image in the elastic band using 
some description of the energetics (ab initio or empir- 
ical force fields). Then for each image, the local tan- 
gent to the path is estimated, and the force defined in 
Eq. H25|l is evaluated for an initial guess of the path. 
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The subsequent minimization for the magnitude of the 
forces with respect to the coordinates of the interme- 
diate images can be carried out with the vefocity Ver- 
let algorithm^^^. Recently, several improvements of the 
original NEB have been proposediS2J^2JS. The NEB 
method has found a wide range of applications in mate- 
rials problems, including cross-slip of screw dislocations 
in metalsi^, diffusion and atomic exchange processes at 
metal and semiconductor surfacesi^^ti^, dissociative ad- 
sorption of molecules on surfacesiS, and contact for- 
mation of metal tips on surfacei^. The NEB method 
has been implemented in many empirical potential and 
ab initio atomistic approachesi^^iiS^. One drawback of 
the NEB method is the difficulty of choosing appropri- 
ate spring constants. A large spring constant requires 
a small time step in the evolution of states, i.e., more 
images along the path. A small spring constant, on the 
other hand, may fail to achieve the desired uniformity 
of the images along the path, and hence may reduce the 
accuracy for the energy barrier. Furthermore, like other 
methods in this category, the NEB method becomes less 
efficient or even inapplicable to systems with very rough 
energy landscapes. 

Realizing the importance of real dynamical pathways. 
Chandler and collaborators have recently proposed meth- 
ods for statistically sampling dynamical paths (MC sam- 
pling of MD trajectories) that do not require the assump- 
tion of TST or the existence of a single, well defined 
transition state or transition patbiSiiSS. In particular, 
no reaction coordinate is needed to study the dynam- 
ics or kinetics of rare transitionsiSi with these methods. 
In a sense, the transition path sampling methods are 
metaphorically akin to throwing ropes over rough moun- 
tain passes in the dark: "throwing ropes" corresponds to 
shooting short real dynamical trajectories; "in the dark" 
implies that the high-dimensional systems are so complex 
that it is generally impossible to visualize the topogra- 
phy of relevant energy surfaces. Although these meth- 
ods are extremely powerful for treating complex systems 
with rough energy landscapes, they are usually computa- 
tionally demanding. In particular, their efficiency usually 
hinges on the ability to produce new accepted paths from 
old ones, thus they have found limited applications so far. 

Recently, an alternative finite temperature string 
method was proposed which represents transition paths 
by their intrinsic parameterization in order to efficiently 
evolve and sample paths in the path spacaiSi. The string 
method performs a constrained sampling of the equilib- 
rium distribution of the system in hyper-planes normal 
to the transition pathways of a coarse-grained potential 
which need not be determined beforehand. The collec- 
tion of the hyper-planes is parametrized by a string which 
is updated sclf-consistently until it approximates locally 
the correct coordinate associated with the reaction event. 
The region in these planes where the equilibrium distri- 
bution is concentrated determines a transition tube in 
configuration space in which a transition takes place with 
high probability. The string method naturally overcomes 



the spring constant problem in the NEB method owing to 
the intrinsic parameterization of the string, and the dis- 
tribution of the replicas along the chain is automatically 
uniform. The method, however, rests on the assump- 
tion that the equilibrium distribution must be localized 
on the iso-surfaces of the reaction coordinate and these 
iso-surfaces can be locally approximated by the hyper- 
planes. If the effective transition tube is highly curved in 
configuration space this approach may have to be modi- 
fied. 

A method for efficiently generating classical trajecto- 
ries with fixed initial and final boundary conditions has 
recently attracted attention because of its conceptual and 
computational simplicitjiiS^. The approach developed by 
Passerone and Parrinello addresses a very general prob- 
lem: given an initial and a final configuration, what 
are the dynamical paths that connect them? Given a 
classical dynamical system described by a set of coordi- 
nates q, and its trajectory q{t) with boundary conditions 
q(0) = qA and q(T) — qs is determined by locating the 
stationary point of the action S: 

S= f Ciqit),qit))dt, (26) 
Jo 

where £ is the Lagrangian C — T — V, and T and V are 
the kinetic and potential energy, respectively. Following 
the work of Gillilan and Wilson the action S can be 
discretized as 




where qo = qA, Qjv — <lB', A = t/N is the time interval 
and the mass is taken as unitarjiiSi. The stationary so- 
lution of this action is the discretized Euler-Lagrangian 
equation and the corresponding trajectory is identical to 
the well-known Verlet trajectory. The novel part of this 
method is to supplement the action with a penalty func- 
tion which favors the energy conserving trajectories: 

e(q„i?)=5 + Ai^(i?, -i?)2. (28) 

Here fi determines the strength with which the energy 
conservation is enforced, Ej is the instantaneous energy 
given by 

i?j = (qj-q,+i)V(2A2) + V(q,) 

and E is its target value. This is motivated by the fact 
that the physical trajectories have to conserve total en- 
ergy. In all the systems studied, it was found that there 
exists a rather large interval of values such that has 
a minimum close to the Verlet trajectories. In order to 
minimize the Q function more efficiently, Passerone et al. 
make the transformation 
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thus automatically satisfying the boundary conditions. 
The advantage of using over is that a,j has a global 
character. In practice, Q is first optimized with respect 
to a relatively small number of a^, thus capturing the 
global features of the trajectory, and then higher fre- 
quency terms are added. Each time a standard conju- 
gate gradient algorithm is used to minimize Q with only 
the evaluation of the forces. The computational scaling 
is therefore linear in the number of degrees of freedom 
rather than quadratic as in other approaches involving 
the Hessian matrix. 

To illustrate the performance of the method, Passerone 
et al. studied a few simple systems. For a one- 
dimensional double well potential, they found that their 
solutions agree very well with the Verlet trajectory but 
without the calculation of the Hessian matrix. The 
second example was a minimization of a trajectory in 
a two-dimensional configurational space, namely in the 
Mueller potential^^^ . Again they found a satisfactory re- 
sult, namely that the trajectories pass exactly through 
the saddle point and the overall behavior of the trajec- 
tories is physical. The last example was to look at a 
process in which the central atom of a seven-atom two- 
dimensional Lennard- Jones cluster migrates to the sur- 
face. In this case, they claim that their calculations re- 
produce the results from the more elaborate method that 
involves transition path sampling procedure^®^. Finally, 
the authors pointed out that their method can be easily 
implemented within the Car-ParrincUo MD approacbiSi, 
offering a powerful tool for the study of problems in chem- 
istry and materials science. The main advantages of this 
method are the fact that it requires only the calculation 
of the forces, its numerical stability and the quality of 
the trajectories. Furthermore, the method lends itself to 
very efficient parallelization, and it can include naturally 
the multiple time-scale approachesiS. 



C. Escaping free-energy minima 

For many systems, the free-energy surface (FES) may 
have multiple local minima separated by large barriers, 
therefore the time-scale that typical MD and MC simula- 
tions can reach is severely limited. Examples of such sys- 
tems include conformational changes in solution, protein 
folding, and many chemical reactions. A large number 
of methods have been developed to overcome the prob- 
lem, some of which were already mentioned in sections 
IIV Al and llV Bl fsee. for instance, the accelerated dynam- 
ics approach"'^'^'', or the dynamical transition path sam- 
pling method—). 

Very recently, a novel and powerful method for explor- 
ing the properties of multi-dimensional FES of complex 
systems has been proposed by Laio and ParrinelloiS. 
This method combines the ideas of coarse-grained dy- 
namics on the FESi22iiii with those of adaptive bias po- 
tential methodsiS*i24. The method allows the system to 
escape from local minima in the FES and at the same 



time achieves a quantitative determination of the FES 
through the integrated process. This method assumes 
that there exist a finite number of relevant collective co- 
ordinates Si,i = l,n where n is a small number, and the 
free energy J^{s) depends on these parameters. The ex- 
ploration of the FES is guided by the generalized forces 
= —d!F/ds\. In order to estimate the forces more 
efficiently, an ensemble of P replicas of the system is 
introduced, each obeying the constraint that the collec- 
tive coordinates have a pre-assigned value Si = s\. The 
coarse-grained dynamics of the collective coordinates is 
defined as follows: 

(30) 

where tr* = s*/Asi and = F^Asi are the scaled col- 
lective coordinates and forces, respectively. Asi is the 
estimated size of the FES well in the direction Si, |^*| is 
the modulus of the n-th dimensional vector and Sa is 
a dimensionless stepping parameter. After the collective 
coordinates are updated using Eq. (|30|l . a new ensemble 
of replicas of the system with values cr*"*"^ is prepared, and 
new forces F^^^ are calculated for the next iteration. At 
the same time the driving forces are evaluated from the 
microscopic Hamiltonian in short standard microscopic 
MD runs. In order to explore the FES more efficiently, 
the forces in Eq. H3U|I arc replaced by a history-dependent 
term: 

t'<t i 

where the height and the width of the Gaussian, W and 
(5cr, are chosen in order to provide a reasonable balance 
between accuracy and efficiency in exploring the FES. 
The component of the forces coming from the Gaussian 
will discourage the system from revisiting the same spot 
and encourage an efficient exploration of the FES. As the 
system diffuses though the FES, the Gaussian potentials 
accumulate and fill the FES well, allowing the system 
to migrate from well to well. After a while the sum of 
the Gaussian terms will almost exactly compensate the 
underlying FES well^^^. 

A typical example of this behavior can be seen in Fig. 
in which a dynamics in Eq. (|30|l is used to explore 
a one-dimensional PES V(s) with three wells. The dy- 
namics starts from a local minimum that is filled by the 
Gaussians in ~ 20 steps. Then the dynamics escapes 
from the well from the lowest energy saddle point, filling 
the second well in ~ 80 steps. The second highest sad- 
dle point is reached in ~ 160 steps, and the full PES is 
filled in a total of ~ 320 steps. Hence, in the case of this 
example, since the form of the potential is known, it can 
be verified that for large t and small (5ct, 

- ^ I^e"^^^ ^ V(s) (32) 

v<t 
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modulo an additive constant. Laio and Parrinello also 
suggest that Eq. holds for a FES, and the free en- 
ergy can be estimated from Eq. 1)32(1 for large The 
efficiency of the method in filling a well in the PES or 
the FES can be estimated by the number of Gaussians 
that are needed to fill the well. This number is propor- 
tional to {\/5a) , where n is the dimensionality of the 
problem. Hence, the efficiency of the method scales ex- 
ponentially with the number of dimensions involved. A 
judicious choice of As^, W and Sa will ensure the right 
balance between accuracy and sampling efficiency, and 
the optimal height and width of the Gaussians can be 
determined by the typical variations of the FES. 

The method has been applied to the study of dissocia- 
tion of a NaCl molecule in water and isomerization of ala- 
nine dipeptide in wateiii^. Overall the method is very ef- 
ficient in exploring the FES of complex systems if the col- 
lective coordinates are chosen judiciously. In particular, 
the topology of a FES can usually be determined by a few 
coarse-grained dynamics steps using "large" Gaussians. 
Subsequently, the qualitative knowledge of the FES can 
be improved using "smaller" Gaussians, effectively re- 
ducing the dimensionality of the problem by exploiting 
the topological information obtained with "large" Gaus- 
sians. As we alluded earlier, the current method assumes 
that the free energy J-{si) depends on a small number of 
collective coordinates Si. However, it is not always obvi- 
ous or possible to identify such collective coordinates for 
complex systems a priori. In the example of isomeriza- 
tion of alanine dipeptide in water, Laio et al. chose the 
dihedral angles $ and as the collective coordinates to 
explore the FES. These authors recognized that the dihe- 
dral angles alone do not provide the complete description 
of the dialanine isomerization reaction, and the real re- 
action coordinates should include the solvent degrees of 
freedom. But their results seemed to reproduce the essen- 
tial features of the FES, therefore the authors concluded 
that the neglected degrees of freedom, although relevant 
for determining the reaction coordinates, are associated 
with small free energy barriers and are sampled efficiently 
during the microscopic dynamics of the dihedral angles 
$ and Despite the success of this particular example, 
identifying a small number of collective coordinates a pri- 
ori., remains challenging within this approach. Moreover, 
the exploration of the FES will be more efficient if an 
adaptive way of determining the parameters of Gaussians 
could be developed. 



D. Other methods 

For systems with a natural disparity in inherent 
time-scales, various multiple-time-step integration algo- 
rithms have been developed to deal with them more 
efficientljiiSiiiSiiZi. One well-known example of such 
strategy is the Born-Oppcnhcimer approximation where 
the electron motion is separated from that of the nuclei 
because of the large disparity between their masses. In 



general, the separation of time-scales occurs when some 
subset of forces present in the system is much stronger 
compared to the rest of the forces while the masses of 
the constituents are about the same. For example, in the 
simulations of the polyatomic liquids with flexible bonds, 
the bond vibrations usually occur at a much faster rate 
than bond translation and rotations. In such systems, 
the configuration space can be divided into fast and slow 
degrees of freedom with the force also being separated 
into fast and slow components. This separation yields a 
set of coupled equations of motion for the evolution of 
the fast and slow degrees of freedom. Instead of solv- 
ing this set of equations simultaneously, multiple-time- 
step integration uses a small time step St to advance the 
fast degrees of freedom n steps holding the slow variables 
fixed. The slow degrees of freedom are then updated us- 
ing a time step nSt. n can be chosen typically between 
5 and 10 in MD simulations of molecules that can be 
described in this way. Furthermore, if an analytic solu- 
tion of high frequency motion can be found, this solution 
can be incorporated into an integration scheme for the 
whole system such that a time step characteristic of the 
slow degrees of freedom can be used and the system can 
be simulated effectively with a much smaller number of 
cyclesi^i. 

Recently, a method based on optimization of the ac- 
tion functional was proposed to extend the time-scale 
of MD simulations by several orders of magnitudeii^. In 
this method, instead of parameterizing the trajectory as a 
function of time, the trajectory is parametrized as a func- 
tion of length. Instead of solving the Newton equations 
in MD simulations, an action term (stochastic difference 
equation with respect to time) is optimized. For acti- 
vated processes the method eliminates the "incubation 
time" , and has proven to be very efficient in the simu- 
lations of biomolecules. It remains to be seen, however, 
if the method can be applied to problems in materials 
science. 



V. CONCLUSIONS 

In this article we have attempted to provide a compre- 
hensive, if not exhaustive, overview of the current status 
of multiscale simulations methods and their applications 
in materials science. We divided the methods that ad- 
dress multiscale problems in the spatial regime into se- 
quential and concurrent. 

The sequential multiscale modeling techniques are in 
general more efficient computationally, but they depend 
on a priori knowledge of physical quantities of inter- 
est, such as the 7-surface in the P-N model, the free- 
energies in the phase-field model, and atomistic local laws 
for mesoscopic DD simulations. The relevance of these 
quantities to the coarse-grained models needs to be care- 
fully examined before the application of the methods. 
Furthermore, these approaches should only be pursued 
when phenomenological theories (such as the P-N model 
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or the phase-field model) are well established, therefore 
the methods are restricted in their range of application. 
In particular, these phenomenological models are often 
associated with the assumption of locality (both in space 
and time). The example of a local approximation in the 
phase-field model is embodied in Eq. (|10|l . which as- 
sumes that part of the energetics of an inhoniogeneous 
system can be written in terms of quantities obtained 
for homogeneous systemsi. Similarly, in the P-N model, 
the 7-energy is assumed to be constant within Aa; dis- 
tance [see Eq. {T))] in order to evaluate the total misfit 
energy. The static approximation (locality in time) for 
dynamical properties is also widely used in phenomeno- 
logical models. The coupling between different scales in 
a sequential approach is usually implicit. A successful 
sequential simulation depends equally on the reliability 
of the phenomenological model and the accuracy of the 
relevant parameters entering the model. 

The concurrent multiscale approaches are much more 
complicated and computationally demanding, but they 
do not require a priori knowledge of physical quantities 
supplied from distinct, lower scale simulations. Further- 
more, concurrent approaches do not depend on any phe- 
nomenological models, therefore they are of more gen- 
eral applicability. Although concurrent approaches are 
more desirable and appealing, the actual problem to be 
attacked must be carefully posed in order to make the 
method practical. The problems that may arise in a con- 
current approach are usually associated with the parti- 
tion of domains in the system. For example, one needs to 
dynamically track the domain boundaries in the MAAD 
simulations and to adapt the EE meshes in the quasi- 
continuum simulations, both of which require additional 
care and computational resources. More importantly, in 
contrast to a sequential method, a "good" hand-shaking 
in a concurrent approach between different domains is 
both challenging and critical. Although some interest- 
ing ideas have been proposed to remedy the problems of 
coupling between different domains, such as the reflec- 
tion of phonons at the domain interfacsi^SiiSSiiSi, there 
is no general consensus on what a proper coupling of do- 
mains is. A general criterion that measures the quality 
of hand-shaking between domains would therefore be de- 
sirable. There is plenty of room for innovative research 
on the issue of domain coupling. General mathematical 
formulations of multiscale problems, including error es- 
timation, may turn out to be very useful for practical 
simulationsiiiiS. In our view, a successful concurrent 
approach usually has to satisfy three conditions: 
(1) Solve the coupled problem (Hamiltonian) accurately 
and efficiently by using ideas such as coarse-graining, 
"bonding energy partition" or "embedding potential" . 



(2) The separate models employed in different domains of 
the system ought to be compatible, i.e., the physical de- 
scription of the system due to the distinct models should 
be as close as possible; 

(3) At each level, the individual model should provide a 
good description of its assigned domain. 

We wish to emphasize the importance of the second con- 
dition which is not in general well recognized. The first 
condition usually guarantees a "smooth" hand-shaking 
across two domains (e.g., the electron density distribu- 
tion or the displacement field varies smoothly across 
the interface), but non-physical charge transfer and/or 
atomic relaxation at the interface could occur if the sec- 
ond condition is not satisfied. Therefore a "smooth" 
hand-shaking does not constitute a "good" hand-shaking, 
and a successful concurrent approach relies on hand- 
shakings that are both mathematically accurate and 
physically consistent. 

In this overview, we have also described a number of 
approaches that strive to extend the temporal scale in 
the modeling and simulation of material properties. We 
categorized these approaches to methods for accelerating 
the dynamics, methods for finding transition paths be- 
tween equilibrium structures, and methods for escaping 
free-energy minima. Although these approaches repre- 
sent very significant developments in the field, the prob- 
lem of linking the time scale of atomic motion and vibra- 
tions (of order femtoseconds) to scales where interesting 
physical phenomena are typically studied (microseconds 
and beyond) is still wide open in many respects. 

Because of the tremendous and continuing progress in 
multiscale strategies, this review is by no means exhaus- 
tive. We hope that we have conveyed the message that 
multiscale modeling is a truly vibrant enterprise of multi- 
disciplinary nature. It combines the skills of physicists, 
materials scientists, chemists, mechanical and chemical 
engineers, applied mathematicians and computer scien- 
tists. The marriage of disciplines and the concomitant 
dissolution of traditional barriers between them, repre- 
sent the true power and embody the great promise of 
multiscale approaches for enhancing our understanding 
of, and our ability to control complex physical phenom- 
ena. 

We acknowledge support from Grant No. E49620-99-1- 
0272 through U.S. Air Force Office for Scientific Research 
(Brown University MURI) and from the Harvard Mate- 
rials Research Science and Engineering Center which is 
funded by the National Science Foundation. We thank 
V.B. Shenoy and A. Laio for providing Fig. |H| and Fig. 
[3 We are grateful to Nick Choly, Weinan E, Paul Mara- 
gakis and Sidney Yip for useful comments and a critical 
reading of the manuscript. 



^ R. Phillips, Crystals, Defects and Micro structures - Mod- 
eling Across Scales (Cambridge University Press, Cam- 
bridge, UK, 2001). 



^ E. Kaxiras and S. Yip (Guest Editors), Curr. Opin. Solid 
State Mater. Sci. 3, 523 (1998) and accompanying articles. 
^ T. Diaz de la Rubia, and V.V. Bulatov (Guest Editors), 



21 



MRS Bull. 26, 169 (2001) and accompanying articles. 
^ S. Yip, Nature Mater., 2, 3 (2003). 

W.M.C. Foulkes, L. Mitas, R.J. Needs, and G. Rajagopal, 

Rev. Mod. Phys. 73, 33 (2001). 
^ A. Szabo, N.S. Ostlund, Modem Quantum Chemistry 

(McGraw-Hill, New York, 1989). 

P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964); 
W. Kohn and L. Sham, Phys. Rev. 140, A1133 (1965). 
* M.C. Payne, M.P. Teter, D.C. Allan, T.A. Arias, and J.D. 

Joannopoulos, Rev. Mod. Phys. 64, 1045 (1992). 
^ S. Goedecker, Rev. Mod. Phys. 71, 1085 (1999). 
1° L. Colombo, Comput. Mater. Sci. 12, 278 (1998). 
" M.Z. Bazant and E. Kaxiras, Phys. Rev. Lett. 77, 4370 
(1996). 

F. Ercolessi and J.B. Adams, Europhys. Lett. 26, 583 

(1994). 

F.F. Abraham, R. Walkup, H. Gao, M. Duchaineau, T. 
de la Rubia, and M. Seager, Proc. Natl. Acad. Sci. USA 
99, 5777 (2002). 

L.P. Kubin and G.R. Canova, Scripta Metall. Mater. 27, 

957 (1992). 

J. P. Hirth, M. Rliee, and H.M. Zbib, J. Comput. Aided 

Mater. Des. 3, 164 (1996). 
^® K.W. Schwarz, Phys. Rev. Lett. 78, 4785 (1997). 
" A. Needleman, Acta Mater. 48, 105 (2000). 
^® V.V. Bulatov, M. Tang, and H.M. Zbib, MRS Bull. 26, 

191 (2001). 

M. Rhee, H.M. Zbib, J.P. Hirth, H. Huang and T. de la 
Rubia, Model. Simul. Mater. Sci. Eng. 6, 467 (1998). 
^° T.J.R. Hughes, The Finite Element Method (Prentice- 
Hall, Englewood Cliffs, NJ, 1987). 

21 P.R. Dawson, E.B. Marin, Adv. Appl. Mech. 34, 77 
(1998). 

22 T. Ohashi, Philo. Mag. Lett. 75, 51 (1997). 

23 D.G. Pettifor, I.I. Oleinik, Phys. Rev. Lett. 84, 4124 
(2000). 

2* D.G. Pettifor, Phys. Rev. Lett. 63, 2480 (1989). 

2® F. Abraham, J. Brougliton, N. Bernstein, and E. Kaxiras, 

Comput. in Phys. 12, 538 (1998). 
2^ J. Broughton, F. Abraham, N. Bernstein, and E. Kaxiras, 

Phys. Rev. B 60, 2391 (1999). 
2^ E. Clementi, Philos. Trans. R. Soc. London, A 326, 445 

(1988). 

2® Computational and Mathematical Models of Microstruc- 
tural Evolution, edited by J.W. BuUard, L.Q. Chen, R.K. 
Kalia, and A.M. Stoneham (Mater. Res. Soc. Proc. 529, 
Warrendale, PA, 1998). 

2® Multiscale Modeling of Materials, edited by V.V. Bula- 
tov, T. Diaz de la Rubia, R. Phillips, E. Kaxiras, and N. 
Ghoniem (Mater. Res. Soc. Proc. 538, Warrendale, PA, 
1999). 

Multiscale Phenomena in Materials - Experiments and 
Modeling, edited by I.M. Robertson, D.H. Lassila, B. 
Devincre, and R. Phillips (Mater. Res. Soc. Proc. 578, 
Warrendale, PA, 2000). 

31 Multiscale Modeling of Materials - 2000, edited by L.P. 
Kubin, J.L. Bassani, K. Cho, H. Gao, and R.L.B. Selinger 
(Mater. Res. Soc. Proc. 653, Warrendale, PA 2001). 

32 G.H. Campbell, S.M. Foiles, H. Huang, D.A. Hughes, 
W.E. King, D.H. Lassila, D.J. Nikkei, T. Diaaz de la 
Rubia, J.Y. Shu, and V.P. Smyshlyaev, Mater. Sci. Eng. 
A251, 1 (1998). 

33 W. E and B. Engquist, Notices of the AMS, 50, 1062 
(2003). 



M.S. Duesbery and G.Y. Richardson, CRC Grit. Rev. 

Soild State Mater. Sci. 17, 1 (1991). 

V. Vitek, Prog. Mater. Sci. 36, 1 (1992). 

B. Joos, Q. Ren, and M.S. Duesbery, Phys. Rev. B 50, 

5890 (1994). 

B. Joos and M.S. Duesbery, Phys. Rev. Lett. 78, 266 

(1997). 

Y. Juan and E. Kaxiras, Philos. Mag. A 74, 1367 (1996). 
J. Hartford, B. von Sydow, G. Wahnstrom, and B.I. 
Lundqvist, Phys. Rev. B 58, 2487 (1998). 
B. von Sydow, J. Hartford, and G. Washnstrom, Comput. 
Mater. Sci. 15, 367 (1999). 

N.I. Medvedeva, O.N. Mryasov, Y.N. Gornostyrev, D.L. 
Novikov and A.J. Freeman, Phys. Rev. B 54, 13506 

(1996) . 

O.N. Mryasov, Y.N. Gornostyrev, and A.J. Freeman, 
Phys. Rev. B 58, 11927 (1998). 

V.V. Bulatov and E. Kaxiras, Phys. Rev. Lett. 78, 4221 

(1997) . 

G. Lu, N. Kioussis, V. V. Bulatov, and E. Kaxiras, Phys. 
Rev. B 62, 3099 (2000); Philos. Mag. Lett. 80, 675 (2000). 
G. Lu, Q. Zhang, N. Kioussis, and E. Kaxiras, Phys. Rev. 
Lett. 87, 095501 (2001). 

G. Lu and E. Kaxiras, Phys. Rev. Lett. 89, 105501 (2002). 
G. Lu, V. V. Bulatov, N. Kioussis, Phys. Rev. B, 66, 
144103 (2002). 

R. Peierls, Proc. Phys. Soc. London 52, 34 (1940). 

F. R.N. Nabarro, Proc. Phys. Soc. London 59, 256 (1947). 
V. Vitek, Philos. Mag. 18, 773 (1968). 

J.D. Eshelby, Philos. Mag. 40, 903 (1949). 

J.P. Hirth and J. Lothe, Theory of Dislocations, 2nd ed. 

(Wiley, New York, 1992). 

S.M. Myers et al. Rev. Mod. Phys. 64, 559 (1992), and 
references therein. 

G. Lu, D. Orlikowski, I. Park, O. Politano, and E. Kaxi- 
ras, Phys. Rev. B 65, 064102 (2002). 

G. Xu, A.S. Argon, and M. Ortiz, Philos. Mag. A 75, 341 

(1997). 

U.V. Waghmare, V. Bulatov, E. Kaxiras, and M.S. Dues- 
bery, Mater. Sci. Eng. A261, 147 (1999). 
U.V. Waghmare, E. Kaxiras, V. Bulatov, and M.S. Dues- 
bery, Model. Simul. Mater. Sci. Eng. 6, 483 (1998). 
P. Peralta, S.A. Maloy, F. Chu, J.J. Petrovic, and T.E. 
Mitchell, Scripta Mater. 37, 1599 (1997). 
J.R. Rice, J. Mech. Phys. Solids 40, 239 (1992); J.R. Rice, 
G.E. Beltz and Y. Sun, in Topics in Fracture and Fatigue, 
ed. A.S. Argon (Springer- Verlag, Berlin 1992). 
Z. Suo, Adv. Appl. Mech. 33, 193 (1997). 
A. Cocks and S. Gill, Adv. Appl. Mech. 36, 81 (1999). 
M. Gurtin, Physica D 92, 178 (1996). 
J. BuUard, E. Garboczi, and W. Carter, J. Appl. Phys. 
83, 4477 (1998). 

D. Fan and L.Q. Chen, Philos. Mag. Lett. 75, 187 (1997). 
L.Q. Chen, C. Wolverton, V. Vaithyanathan, and Z.K. 
Liu, MRS Bull. 26, 197 (2001). 

V. Vaithyanathan, C. Wolverton, and L.Q. Chen, Phys. 
Rev. Lett. 88, 125503 (2002). 

Y.Z. Wang, L.Q. Chen, and A.G. Khachaturyan, 

in Computer Simulation in Materials Science 

Nano/Meso/Macroscoptc Space and Time Scales, 

NATO ASI Series, edited by H. Kirchner, L. Kubin, and 

V. Pontikis (He d'Oleron, Prance, 1995). 

L.Q. Chen and Y.Z. Wang, JOM 48, 13 (1996). 

D.Y. Li and L.Q. Chen, Acta mater. 46, 2573 (1998). 



22 



™ A.G. Khachaturyan, Theory of Structural Transforma- 
tions in Solids (John Wiley, Now York, 1983). 
J.W. Cahn, Acta Mctall. 9, 795 (1961). 

'^^ P.O. Hohenberg and B.I. Halperin, Rev. Mod. Phys. 49, 
435 (1977). 

''^ A. Zunger, in Statics and Dynamics of Alloy Phase Trans- 
formations, NATO ASI Series, edited by P.E.A. Turchi 
and A. Gonis (Plenum Press, New York, 1994) p. 361. 

""^ D. de Fontaine, Solid State Phys. 47, 33 (1994). 

''^ C. Wolverton, Philos. Mag. Lett. 79, 683 (1999); Model. 
Simul. Mater. Sci. Eng. 8, 323 (2000). 

''^ D.B. Laks, L.G. Ferreira, S. Froyen, and A. Zunger, Phys. 
Rev. B 46, 12587 (1992). 

''^ L. Kaufman and H. Bernstein, Computer Calculation of 
Phase Diagrams with Special Reference to Refractory Met- 
als (Academic Press, New York, 1970). 
N. Saunders and A.P. Miodownik, CALPHAD: A Com- 
prehensive Guide (Pergamon Press, Oxford, 1998). 

'■^ P.J. Spencer, MRS Bull. 24, 18 (1999). 

®° D. Kandel and E. Kaxiras, Solid State Phys. 54, 219 

(2000) . 

P. Kratzer and M. Scheffler, Compt. Sci. Eng. 3, 16 

(2001) . 

W. Cai, V.V. Bulatov, J.F. Justo, S. Yip, and A.S. Argon, 
Phys. Rev. Lett. 84, 3346 (2000). 

C. Ratsch, M.F. Gyure, R.E. Caflisch, M. Petersen, M. 
Kang, J. Garcia, and D.D. Vvedensky, Phys. Rev. B 65, 
195403 (2002). 

®* M.F. Gyure, C. Ratsch, B. Merriman, R.E. Caflisch, S. 
Osher, J. Zinck, and D.D. Vvedensky, Phys. Rev. E 58, 
6927 (1998). 

*5 S. Oshcr and J. A. Sothian, J. Compt. Phys. 79, 12 (1988). 
C.W. Shu and S. Osher, J. Compt. Phys. 83, 32 (1989). 
C.S. Deo, D.J. Srolovitz, W. Cai and V.V. Bulatov, to be 
published. 

V. Bulatov, F.F. Abraham, L. Kubin, B. Devincre and S. 
Yip, Nature 391, 669 (1998). 

A. Needleman and E. Van der Giessen, MRS Bull. 26, 211 

(2001) . 

®° F. Stillinger and T.A. Weber, Phys. Rev. B 31, 5262 
(1985). 

N. Bernstein and E. Kaxiras, Phys. Rev. B 56, 10488 
(1997). 

"2 K. Ohsawa and E. Kuramoto, J. Appl. Phys. 86, 179 

(1999) . 

"^^ F.F. Abraham and H. Gao, Phys. Rev. Lett. 84, 3113 

(2000) . 

S.J. Zhou, P.S. Lomdahl, R. Thomson, and B.L. Holian, 
Phys. Rev. Lett. 76, 2318 (1996). 

B. L. Holian and R. Ravelo, Phys. Rev. B 51, 11275 
(1995). 

P. Gumbsch, S.J. Zhou, and B.L. Hohan, Phys. Rev. B 

55, 3445 (1997). 

S.J. Carroll, P.D. Nellist, R.E. Palmer, S. Hobday, and R. 

Smith, Phys. Rev. Lett. 84, 2654 (2000). 

M. Moseler, J. Nordiek, and H. Haberland, Phys. Rev. B 

56, 15439 (1997). 

^3 N. Choly and E. Kaxiras, SoUd State Comm. 121, 281 

(2002) . 

S. Watson, and E. Carter, Comput. Phys. Comm. 128, 
67 (2000). 

S. KohlhofF, P. Gumbsch, H.F. Fischmeister, Philos. Mag. 

A 64, 851 (1991). 

P. Gumbsch, J. Mater. Res. 10, 2897 (1995). 



E.B. Tadrnor, M. Ortiz, and R. Phillips, Philos. Mag. A 
73, 1529 (1996). 
1°"* E.B. Tadrnor, R. Phillips, and M. Ortiz, Langmuir 12, 
4529 (1996). 

V.B. Shenoy, R. MiUer, E.B. Tadrnor, R. PhiUips, and M. 

Ortiz, Phys. Rev. Lett. 80, 742 (1998). 

R. Miller, E.B. Tadmor, R. Phillips, and M. Ortiz, Model. 

Simul. Mater. Sci. Eng. 6, 607 (1998). 
1°'' R. Miller, M. Ortiz, R. Phillips, V. Shenoy, and E. B. 
Tadmor, Eng. Fracture Mech. 61, 427 (1998). 

D. Rodney and R. Phillips, Phys. Rev. Lett. 82, 1704 
(1999). 

V.B. Shenoy, R. Miller, E.B. Tadmor, D. Rodney, R. 
Phillips, and M. Ortiz, J. Mech. Phys. Solid 47, 611 

(1999) . 

"° E.B. Tadmor, R. Miller, R. Phillips, and M. Ortiz, J. 
Mater. Res. 14, 2233 (1999). 

V.B. Shenoy, R. Phillips, and E.B. Tadmor, J. Mech. 
Phys. Solid 48, 649 (2000). 

G.S. Smith, E.B. Tadmor, and E. Kaxiras, Phys. Rev. 
Lett. 84, 1260 (2000). 
"3 J. Knap and M. Ortiz, J. Mech. Phys. Solid 49, 1899 

(2001) . 

G.S. Smith, E.B. Tadmor, N. Bernstein, and E. Kaxiras, 
Acta Mater. 49, 4089 (2001). 

V. Shenoy, V. Shenoy, and R. Philhps, in Multiscale Mod- 
eling of Materials, edited by V.V. Bulatov, T. Diaz de la 
Rubia, R. Phillips, E. Kaxiras, and N. Ghoniem (Mater. 
Res. Soc. Symp. Proc. 538, Warrendale, PA, 1999), p. 
465. 

A.H. King and D.A. Smith, Acta Crystallogr. A 36, 335 

(1980). 

E. B. Tadmor, G.S. Smith, N. Bernstein, and E. Kaxiras, 
Phys. Rev. B 59, 235 (1999). 

"8 E.B. Tadmor, U.V. Waghmare, G.S. Smith, and E. Kaxi- 
ras, Acta Mater. 50, 2989 (2002). 

"® R.E. Rudd and J.Q. Broughton, Phys. Rev. B 58, R5893 
(1998). 

^™ R.E. Rudd and J.Q. Broughton, Phys. Stat. Sol. 217, 251 

(2000) . 

R.E. Rudd and J.Q. Broughton (unpublished). 

S. Curtarolo and G. Ceder. Phys. Rev. Lett. 88, 255504 

(2002) . 

A. A. Migdal, Zh. Eksp. Tcor. Fiz. 69, 1457 (1975); L.P. 
Kadanoff, Ann. Phys. (N.Y.) 100, 259 (1968). 

L.E. Shilkrot, R.E. Miller and W.A. Curtin, Phys. Rev. 
Lett. 89, 025501 (2002). 

E. van der Giessen and A. Needleman, Model. Simul. 
Mater. Sci. Eng. 3, 689 (1995). 

N. Bernstein and D. Hess, Mat. Res. Soc. Symp. Proc. 
653, Z2.7.1 (2001). 
^'^'^ E. Lidorikis, M. Bachlechner, R.K. Kalia, A. Nakano, 
P. Vashishta, and G.Z. Voyiadjis, Phys. Rev. Lett. 87, 
086104 (2001). 

C. Woodward and S.I. Rao, Phys. Rev. Lett. 88, 216402 
(2002). 

J. Gao, Rev. Compt. Chem. 7, 119 (1996). 

Q. Cui, M. Elstner, E. Kaxiras, T. Frauenheim, and M. 

Karplus, J. Phys. Chem. B 105, 569 (2001). 

M. Elstner, D. Porezag, G. Jungnickel, J. Eisner, M. 

Haugk, T. Frauenheim, S. Suhai, and D. Seifert, Phys. 

Rev. B 58, 7260 (1998). 

B. R. Brooks, R.E. Burccoleri, B.D. Olafson, D.J. States, 
S. Swaminathan, and M. Karplus, J. Compt. Chem. 4, 



23 



187 (1983). 

N. Govind, Y.A. Wang, A.J.R. da Silva, and E.A. Carter, 
Chem. Phys. Lett., 295, 129 (1998). 
^^"^ N. Govind, Y.A. Wang, and E.A. Carter, J. Chem. Phys. 
110, 7677 (1999). 

T.A. Wesolowski and A. Warshel, J. Phys. Chem. 97, 8050 

(1993) . 

T. Kluner, N. Govind, Y.A. Wang, and E.A. Carter, Phys. 
Rev. Lett. 86, 5694 (2001). 

A.F. Voter, F. Montalenti, and T.C. Germarm, Anrm. 

Rev. Mater. Res. 32, 321 (2002). 

A.F. Voter, J. Chem. Phys. 106, 4665 (1997). 

A. F. Voter, Phys. Rev. Lett. 78, 3908 (1997). 
R. Marcelin, Ann. Phys. 3, 120 (1915). 

"1 E. Wigner, Z. Phys. Chem. B 19, 203 (1932). 
H. Eyring, J. Chem. Phys. 3, 107 (1935). 
P. Hanggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. 
62, 251 (1990). 

J.B. Anderson, Adv. Chem. Phys. 91, 381 (1995). 
"5 G.H. Vineyard, J. Phys. Chem. SoUd 3, 121 (1957). 
"'5 A.F. Voter, Phys. Rev. B 57, R13985 (1998). 
^^"^ M. S0rensen and A.F. Voter, J. Chem. Phys. 112, 9599 

(2000). 

R. Elber and M. Karplus, Chem. Phys. Lett. 139, 375 
(1987). 

^■^^ H. Jonsson, G. Mills, and K.W. Jacobsen, in Computer 
Simulation of Rare Events and Dynamics of Classical 
and Quantum Condensed-Phase Systems, edited by B.J. 
Berne, G. Ciccotti, and D. Coker (World Scientific, Sin- 
gapore, 1998). 

G. Henkelman, G. Johannesson, and H. Jonsson, in- 
Progress on Theoretical Chemistry and Physics, edited 
by S. D. Schwartz ( Kluwar Academic Publishers, 2000). 

H. C. Anderson, J. Chem. Phys. 72, 2384 (1980). 

G. Henkelman, B. Uberuaga and H. Jonsson, J. Chem. 

Phys. 113, 9901 (2000). 
153 p Maragakis, S. Andreev, Y. Brumer, D. Reichman, and 

E. Kaxiras, J. Chem. Phys. 117, 4651 (2002). 
^^"^ T. Rasmussen, K.W. Jacobsen, T. Leffers, O.B. Petersen, 

S.G. Srinivasan, and H. Jonsson, Phys. Rev. Lett. 79, 

3676 (1997). 

M. Villarba and H. Jonsson, Surf. Sci., 317, 35 (1994). 

B. P. Uberuaga, M. Leskovar, A. P. Smith, H. Jonsson, 
and M. Olmstead, Phys. Rev. Lett., 84, 2441 (2000). 

G. Mills and H. Jonsson, Phys. Rev. Lett., 72, 1124 

(1994) . 

M.R. Sorensen, K.W. Jacobsen, and H. Josson, Phys. Rev. 



Lett., 77, 5067 (1996). 

C. Dellago, P.G. Bolhuis, F.S. Csajka, and D. Chandler, 
J. Chem. Phys. 108, 1964 (1998). 

P.G. Bolhuis, C. Dellago, and D. Chandler, Faraday Dis- 
cuss. Chem. Soc. 110, 421 (1998). 

P.G. Bolhuis, D. Chandler, C. DeUago, and P.L. Geissler, 

Annu. Rev. Phys. Chem. 53, 291 (2002). 

W. E, W. Ren and E. Vanden-Eijnden, to be pubUshed. 

D. Passerone and M. Parrinello, Phys. Rev. Lett. 87, 
108302 (2001). 

R.E. Gillilan and K.R. Wilson, J. Chem. Phys. 97, 1757 
(1992). 

^'^^ K. Mueller, Angew. Chem. 19, 1 (1980). 

C. Dellago, P.G. Bolhuis, and D. Chandler, J. Chem. 
Phys. 108, 9326 (1998). 

R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 

(1985) . 

V. Zaloj and R. Elber, Comput. Phys. Commun. 128, 118 

(2000) . 

^'^^ A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. USA 

99, 12562 (2002). 
'^^^ l.G. Kevrekidis, C. Theodoropoulos, and C.W. Gear, 

Computers and Chemical Engineering (to be published). 

W. E and B. Engquist, Comm. Math Sci. 1, 87 (2003). 

W. E, B. Engquist, and Z. Huang, Phys. Rev. B 67, 

092101 (2003). 

T. Huber, A.E. Torda, and W.F. van Gunsteren, J. Com- 
put. Aided Mol. Design 8, 695 (1994). 
F. Wang and D.P. Landau, Phys. Rev. Lett. 86, 2050 

(2001) . 

R.D. SwindoU and J.M. Haile, J. Comput. Phys. 53, 289 
(1984). 

^''^ O. Teleman and B. Jonsson, J. Comput. Chem. 7, 58 

(1986) . 

^'^'^ M.E. Tuckerman, G.J. Martyna, and B.J. Berne, J. Chem. 

Phys. 93, 1287 (1990). 
'^^^ R. Elber, A. Ghosh, and A. Cardenas, Acc. Chem. Res. 

35, 396 (2002). 

W. Cai, M. de Koning, V.V. Bulatov, and S. Yip, Phys. 
Rev. Lett. 85, 3213 (2000). 

W. E and Z. Huang, Phys. Rev. Lett. 87, 135501 (2001). 
W. E and Z. Huang, J. Comput. Phys. 182, 234 (2002). 
A web site with useful information related to 
the quasicontinuum method can be found at 
http://www.qcmethod.com where the quasicontin- 
uum codes are also available to download. 



24 




FIG. 1: A schematic illustration of spacial and temporal scales achievable by various simulation approaches. The scales are in 
centimeters for the length dimension and seconds for the time dimension, both logarithmic. 
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FIG. 2: A schematic illustration showing an edge dislocation in a lattice. The partition of the dislocated lattice into linear 
elastic region and nonlinear atomistic region allows a multiscale treatment of the problem. 
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FIG. 3: The 7-surface (J/ra^) for displacements along a (111) plane for (a) pure Al and (b) Al+H systems. The corners of the 
plane and its center correspond to identical equilibrium configurations, i.e., the ideal lattice. The two surfaces are displayed in 
exactly the same perspective and on the same energy scale to facilitate comparison of important features. 
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FIG. 4: Illustration of the three different levels of simulation in the level-set multiscale approach of surface growth, (a) The 
macroscopic scale, in which island borders are continuous lines separating heights at different levels (the levels along a particular 
cross section are shown schematically, labeled 0, 1 and 2); this scale is treated with the level-set method, (b) The mesoscopic 
scale, where the features of the island edges contain some information about the underlying atomic lattice, indicated here 
as the small straight lines that define step directions consistent with atomic positions; this scale is treated with the kinetic 
Monte Carlo approach, (c) The microscopic scale, where the individual degrees of freedom are explicitly included. The step 
is determined by the positions of atoms in two terraces, the ones on the upper terrace shown as white larger circles while the 
ones on the lower terrace shown as shaded smaller circles; this scale is treated by atomistic (possibly ab initio) methods. All 
views in this schematic representation are top views (see text for details). 
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FIG. 5: Geometrical decomposition of a Si slab with a small crack into different dynamic regions in the a MAAD simulation: 
(a) The system at the macroscopic scale, which is modeled as a continuum using finite elements (FE), except for the region near 
the crack outlined in dashed line, (b) The mesoscopic scale, treated atomistically with interatomic potentials and molecular 
dynamics (MD) with the atoms indicated by green circles, except for the region in the immediate neighborhood of the crack, 
outlined in dashed line, (c) The microcopic scale, treated atomistically with forces derived from quantum mechanical calculations 
with the atoms indicated by red circles. The hand shaking regions between FE and MD and between MD and TB are also 
shown schematically (see text for details). 
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FIG. 6: Example of a rnultiscalc simulation using the quasicontinuum method, (a) Finite-element mesh used to model 
dislocation-GB interaction. The surface marked AB is rigidly indented in order to generate dislocations at A (distance in 
A), (b) Snapshots of atomic positions at different stages in the deformation history. Absorption of the first pair of dislocations 
at the GB results in a step, while the second pair form a pileup. 




FIG. 7: Time evolution of the sum of a one-dimensional model potential V(a) and the accumulating Gaussian terms of Eq. 
1311 1. The dynamic evolution (thin lines) is labelled by the number of dynamical iterations in Eq. 1301 . The starting potential 
(thick line) has three minima and the dynamics is initiated in the second minimum. 



